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ABSTRACT 


The problem of ameliorating the discomfort of passengers on a 
large air transport subject to flight disturbances is examined. The 
longitudinal dynamics of the aircraft, including effects of body 
flexing, are developed in terms of linear, constant coefficient 
differential equations in state variables. A cost functional, 
penalizing the rigid body displacements and flexure accelerations 
over the surface of the aircraft is formulated as a quadratic form. 
The resulting control problem, to minimize the cost subject to the 
state equation constraints, is of a class whose solutions are well 
known. The feedback gains for the optimal controller are calculated 
digitally, and the resulting autopilot is simulated on an analog 
computer and its performance evaluated. 
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INTRODUCTION 


The design of an autopilot system to improve the ride quality of 
transport aircraft is a formidable task. The system commonly accepts 
data from a number of sensors placed at various points throughout 
the airframe, and outputs command signals to the aircraft control sur- 
faces. The autopilot designer must synthesize this multi-input/output 
system and evaluate the performance of the design in terms of the ride 
quality at all seats of the airplane. Conventional methods of control 
system design (Bode, Root Locus, etc.) readily handle questions on 
stability and transient response, but are generally insufficient to 
provide the subtle insights necessary to determine how gain ratios or 
feedback paths should be altered to improve the overall ride quality. 

In the work presented here, it will be shown that this autopilot 
design problem can be formulated realistically as an optimal control 
problem. A solution will be demonstrated for a typical supersonic 
transport aircraft and the resulting design evaluated by simulation 
on an analog computer. 

The first attempt at development of a ride quality autopilot was 
the result of flight test experience with the XB-70 experimental bomber, 
where high sensitivity to input gusts in the forward fuselage area 
resulted in an uncomfortable acceleration environment for the crew. 

(This problem has since been lessened by the implementation of the 
ILAF control system on the XB-70) . While this situation remains 
undesirable in a military aircraft, it is intolerable in a civilian 
transport, where the concept of ride quality becomes exceedingly 
important. 

Aircraft flexing is only one detractor from passenger ride comfort. 
The characteristic vertical motion of the airplane , resulting from 
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excitation of the short period mode, adds to rider discomfort. 

While the total displacements due to airframe flexure may be small, 
typically the accelerations due to the flexing occur at frequencies 
which are most annoying to human subjects. On the other hand, the 
large vertical motion of the aircraft, while not contributing 
significantly to the disturbing accelerations, may cause dizziness 
or motion sickness. Any attempt to control the ride quality of an 
airplane must consider these two differing contributors. 

The Boeing swept-wing SST Transport prototype (B-2707) was 
chosen as a model for this investigation. As a proposed civilian 
airliner, ride comfort is important, and the design is such to in- 
duce greater fuselage flexibility than in previous commercial air- 
planes, and even the XB-70. It should be pointed out, however, that 
unlike the XB-70, the B-2707 flexure modes are well behaved in that 
flexure disturbances are not concentrated in either forward or aft 
fuselage sections. 

As the attempt of this thesis is to investigate the control 
problem associated with minimizing ride discomfort of a large flexible 
aircraft, rather than to design a particular control system for a 
particular airplane, an idealized and simplified mathematical model 
of the B-2707 was chosen. The dimensions of the model are almost 
identical to those of the B-2707, however, the mass distribution of 
the model was assumed constant over the fuselage and constant over the 
lifting surfaces. This two phase distribution allowed for the mass 
and center of gravity of the model to be consistent with those of 
the B-2707. The effective deviation from the actual aircraft is most 
probably small (whenever the mass distribution arises, it is smoothed 
via integration) , yet the freedom from following minute detail of 
the aircraft allows for more insight into the structure of the problem 
at hand. The idealized model is shown in Figure 1. Numerical values 
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relating to physical aspects of the model were either given in the 
Boeing 2707 Airframe Design Report (Ref. 2) , or else were derived to 
conform to some expected behavioral characteristic of the B-2707. In 
either case, the origin of numerical data will be made clear as they 
arise in the text. 

Before one can tackle the problem of flexure control with the 
expectation of manageable results, vast simplifications must be 
achieved. The problem is of a distributed nature; the vibrating 
structure has an infinity of natural modes. The complex dynamic 
interactions of subassemblies, wings, fuselage, engine nacelles, as 
well as concomitant forms of flexing, tortional vs. longitudinal, compli- 
cate the problem to extraordinary degree. Evidently, it would be 
advantageous to effect both a decoupling between flexing of sub- 
assemblies and an approximation of the distributed system by a finite 
number of modes. 

In the succeeding chapter, a model for the flexure dynamics of 
the idealized aircraft will be developed. The particular model will 
obviously depend critically upon the elastic and inertial properties 
of the structure. The flexure analysis in this paper ignores the 
effect of aerodynamics upon the elastic properties of the vehicle. 

The result is a model in which the airplane vibrates according to a 
set of uncoupled modes. If the local changes in aerodynamic loading 
due to aircraft flexure were included, then the forcing terms, for 
example, would no longer be independent of the motion of the aircraft, 
and the resulting equations of vibration would become inexorably 
coupled. However, as the displacements of the elastic structure due 
to flexing are small, the deviations from nominal aerodynamic forces 
will be correspondingly small. 
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E = elevator 


Figure 1. Idealized SST Model (all dimensions in inches) 






It should be understood that large differences in the flight 
condition may change the nominal aerodynamic loading on the aircraft 
structures greatly. Thus, the model derived in this investigation, 
based on conditions of constant speed, constant altitude, steady 
state flight may be inappropriate for other flight conditions, such 
as steep climbing. 

Control surfaces available to the controller will consist solely 
of the elevators.* In the swept wing supersonic phase of flight, 
the elevators represent the major normal force producing control sur- 
faces. They are generally effective in producing pitching moments 
and in tailess aircraft, (In the sweptback wing configuration, the 
B-2707 is essentially a tailess aircraft.) adept also in producing 
left. In addition, the position of the elevators does not coincide 
with any nodes of the flexing motion, hence forces applied through 
the elevators will influence all the deflection modes. 


♦Although only the elevators are considered in this development, 
other control effectors (flaps, spoiler, etc.) could be readily 
included. 
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CHAPTER 1 

T he Physical Problem 

In this chapter, a model of the airplane will be developed suitable 
for a control formulation of the problem. Section A will explore the 
flexure dynamics of the idealized aircraft. The problem of rigid body 
motion and the aerodynamic properties involved are discussed in 
Section B. Control elements and control forces are examined in Section 
C and a suitable cost functional is developed in Section D. In addi- 
tion, supplemental information and detailed calculations have been 
inserted as appendices and are referenced throughout the chapter. 

A. Flexure Dynamics 
A. 1. Introduction 

Although the modal response of an elastic airframe is a distributed 
parameter system, intuitive considerations demand that the response at 
higher frequencies become vanishingly small lest a bounded energy 
input give rise to an unbounded energy output. An in depth analysis 
of the B-2707 confirms this and suggests that only the first five 
vibration modes can be considered at all significant. (Ref. 2, p.52) 
Further data shows that although the frequency response of the airplane 
to inputs at the fifth modal frequency (5.09 cps) is not insignificant, 
the harmonic content of realistic atmosphere gust models is appreciably 
lower. Hence in flight conditions subject to random gust disturbances, 
the excitation of the fifth mode is minimal, and it suffices to examine 
just the first four modes. The natural vibration characteristics 
corresponding to these modes are shown in Figure 2. and refer to 
the lightweight (375,000 lb.) configuration of the airplane. 

It is instructive to examine closely the natural modes, for 
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their particular shapes shed much light on the modeling problem. 

From Figure 2, certain modal characteristics stand out. The first 
and fourth modes (1.47 and 3.81 cps, respectively) consist mainly of 
fuselage flexing and tend to be sinusoidal in shape. The wing and 
tail act to limit the motion of the aft section of the fuselage in 
the first mode, as these lifting surfaces vibrate in tandem. In the 
fourth mode however, the lifting surface motion is seen to be of a 
scissoring nature, and while still introducing damping into the vi- 
brating system, has a minimal effect on the range of motion exhibited 
by the rear sections of the passenger compartment. The second 
mode (2.25 cps) consists almost entirely of lifting surface motion. 
Here again, the wing and tail scissor, inducing scarcely any motion 
along the fuselage. Mode three (2.56 cps) is again primarily lifting 
surface flexure, but in this case, wing and tail motion is in tandem, 
and the effect is greater fuselage flexing than in the former case. 

Judging from Figure 2 and the preceeding discussion, the decompo- 
sition of the airframe into independent flexing subassemblies would 
seem at first glance a profitable avenue of attack. Through a proper 
decomposition, it should be possible to model the subassemblies with 
a degree of homogeneity lacking in the entire aircraft. For example, 
we might want to divorce the wings from the fuselage in our model, 
considering each as an independent structure, and superpositioning 
solutions with appropriate boundary conditions. Whatever the decompo- 
sition, we would demand that the autonomously derived solutions ade- 
quately conform to the expected mode shapes in Figure 2. 

Therein is where the problem lies. The boundary conditions which 
need to be matched up when the components are reassembled are time 
varying. Whatever simplifications occur by the decomposition are more 
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Figure 2. Natural Vibration Characteristics, B-2707 
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than counterbalanced by the complications of the boundary conditions. 
Yet the idea is not a totally fruitless endeavor. The investigation 
of the flexure of the independent subassemblies, if not yielding 
closed form solutions, nevertheless sheds much light on the dominant 
motions of the flexing aircraft. Even more important, the theory of 
long, slender elastic beams which is often brought into play to model 
the various airframe components, extends quite naturally to problems 
of elastic plates. It is the theory of elastic plates which will allow 
us a framework for a simple, yet highlv applicable model for airframe 
flexure . 

For the reasons outlined above, we will take a small digression 
and examine briefly the flexure dynamics of the dominant subassembly, 
the fuselage. Any model of the fuselage dynamics must satisfactorily 
account for the first and fourth flexing modes of Figure 2, as these 
are seen to be composed primarily of fuselage vibrations. A seemingly 
reasonable and yet simple choice would be to equate the fuselage to a 
long slender beam of appropriate mass distribution and stiffness. 
Disregarding the extreme forward and aft fuselage sections, there 
is little reason to expect anything other than a fairly uniform mass 
distribution along the length of the fuselage. This model yields 
mode shapes (Figure 3) which are very similar to the manufacturer's 
vibration characteristics, modes one and four. In addition, the beam 
theory predicts a ratio of 2.77 between the first two normal mode 
frequencies. The manufacturer's data on the B-2707 has shown the two 
fuselage dominant modes to be separated by a factor of 2.56, in close 
accord with the uniform beam theory. 

The application of the beam theory provides us with a model which 
is well understood, and soluble. As the rotary inertia and transverse 
shear deformations of the airplane fuselage are negligible in the 
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assumed flight configuration, these terms may be omitted in the beam 
model, yielding a solution which can be exhibited in closed form. 

Appendix A, Sections 1 and 2, contain derivations of both the differential 
and integral equations of motion for the slender beam. Both analyses 
are included as they represent different mathematical approaches 
to the same problem, and tend to complement each other. The reader 
should note the inherent similarities in the derivations of the inte- 
gral equations of motion for the simple beam and the unrestrained 
aircraft. The latter is given in Appendix 1, Section 3. 

A. 2 Free Vibration of the Unrestrained Airplane 

The two dimensional theory derives its impetus from the simple 
beam analysis. It assumes that a three dimensional aircraft is com- 
pressed into a flat elastic plate, rigid in its own plane and able 
to execute small displacements in vertical position, pitch and roll in 
addition to elastic motion. Although these assumptions do involve a 
loss of generality, they allow us to draw upon previously established 
beam theory in developing a simple yet adequate model for aircraft 
aeroelasticity . 

Consider an airplane idealized as an elastic plate and a super- 
imposed coordinate system, so that the aircraft is in the x-y plane, 
with the x-axis along the length of the fuselage. Furthermore, the 
coordinate axes coincide with the principal inertial axes of the air- 
plane and the origin of the system is at the center of gravity of the 
vehicle. As the structure is vibrating freely, the coordinate axes 
may be considered fixed in space. 

Let us denote the displacement of the elastic surface along the 
z-axis by w(x,y,t). Then in the absence of external forces, moment 


18 



equilibrium about the x and y axes together with force equilibrium 
demand that : 


//w(x,y,t) p (x,y)dxdy = 0 
JJ w(x,y,t)xp (x,y)dxdy = 0 
fly (x,y ,t)yp (x,y)dxdy = 0 


( 1 - 1 ) 


( 1 - 2 ) 


(1-3) 


where p(x,y) represents the mass distribution over the surface of the 
airplane. Additionally, we obtain an equation relating the elastic 
and inertial forces 


w(x,y,t) - w(0, 0 ,t) 


3w(0,0,t) _ „ 3 w(0,0,t) 
3x y 5y 


(1-4) 



C(x,y?4,n)p(5,u)w(5,u,t)d£dn 


where C(x,y;?,ri) is a two dimensional influence function, measuring the 
deflection at a point <x,y) due to a unit force applied at the point 
(£,ri) when the origin of the x-y plane is clamped. Equation (1-4) 
is well motivated by the development of the integral equations of 
motion for a simple beam in Appendix A, Section 2. 

The equations of motion for the elastic plate are separable, 
i.e. they admit a solution of the form 


w(x,y,t) = W(x,y)T(t) 


(1-5) 


Upon the introduction of (1-5) and some simplification, the 
equations can be rewritten in the form: 
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T + w 2 T = 0 


( 1 - 6 ) 


W (x,y) 


ui 



x,y;C,n)W(C,n)p (£ ,n)d£dn 


(1-7) 


where 

G(x,y;£,n) =C(x,y;£,g) ~ffc( r,s;£,q) [^ + ^ + ^-] p (r ,s)drds . (1-8) 

S y x 

Equation (1-8) defines the influence function of the unrestrained air- 
craft. As in the case of the unrestrained beam, the integral term in 
Eq. (1-8) has the effect of subtracting off the rigid body motion of 
the aircraft. Equations (1-6) and (1-7) can be solved for a family 
of solutions, each pair T\ and VT corresponding to a particular value 
of a). The functions satisfying Eq. (1-7) are the natural mode shapes 
of the system and represent the characteristic shapes of the distor- 
tions of the elastic surface. The derivation of the equations of 
free vibration are given in Appendix A, Section 3. In addition, the 
mode shapes are orthogonal to each other with respect to the mass 
distribution function p(x,y): 


/£w m W n Pdxdy = 0 


(1-9) 


This most important characteristic is proved in Appendix A, Section 4. 


A. 3 Forced Motion of the Unrestrained Aircraft 

Before the problem of forced motion of the aircraft can be tackled, 
it is necessary to redefine the meaning of the displacement function 
w(x,y,t). For the situation of free vibration, w(x,y,t), represents 
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the displacement with respect to the airplane principal inertial 
axes, which may be assumed stationary in space. We are afforded the 
opportunity to measure the elastic motion alone by subtracting out 
the gross rigid body motion. Reference to Equation (1-7) will indi- 
cate that the rigid body mode shapes for the freely vibrating beam are 
indeed vacuous. For forced motion, it is no longer possible to pre- 
calculate and subtract the rigid body displacements, hence w(x,y,t) 
must represent the total displacement of the elastic surface including 
gross translations and rotations. 

As is customary when considering only longitudinal dynamics, we 
restrict our study to an elastic airframe free to pitch and plunge, 
but restrained against rolling. F z (x,y,t) will represent the applied 
force in the z-direction per unit surface area. By simple equilibrium 
of forces and moments we have: 


//« (x,y ,t) p (x,y)dxdy = ffo z (x,y,t)dxdy 


( 1 - 10 ) 


ffiHx, y,t) p (x,y)xdxdy = Jf F z (x,y,t)xdxdy . (1-11) 


Analogous to Equation (1-4) the relation between inertial and elastic 
forces becomes 


w(x,y,t) - w (0 , 0 , t) - x 3W ^°' 


= //c(x,y;r,snF z <r,s,t) - w (r , s , t) p (r , s) ] drds . (1-12) 


As in the case of the simple beam [See Appendix A, Section 1 or 2 ] 
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we represent the solution to the equations above as an infinite sum 
of normalized mode shapes and normal (or generalized) coordinates. 


CO 

w (x , y , t) = l <f> (x,y)E (t) 
n= 1 


(1-13) 


The normalized mode shapes <f)j (x,y) are simply the natural mode 
shapes given by the solution of Equation (1-7) multiplied by a parti- 
cular constant. Different procedures for normalizing the mode shapes 
are in wide use. The manufacturer's mode shapes in Figure 2 have 
been scaled so that the maximum value is unity. Another method em- 
ployed is 


M, 


ffs'i*** 


1/2 


W. 

3 


(1-14) 


The use of M. , the aircraft total mass, above is arbitrary as any 
characteristic of the aircraft may be used. The result of normaliza- 
tion by Eq. (1-14) is that the integral over the surface of the squared 
modal deflection weighted by the mass distribution is a constant. This 
integral, the generalized mass, arises in the resultant equations of 
motion. The actual normalization technique used is arbitrary, and in 
fact every mode can have a different and independently derived multi- 
plicative constant without effecting the resultant calculated motion. 

As the displacement function w(x,y,t) consists of the total 
motion of the elastic airframe, one normal mode must illustrate the 
translation of the center of gravity and another the pitching rotation 
about the c.g. Thus we designate 


<)> 1 (x,y) =1; oj 1 = 0 


(1-15) 


<J> 2 (x,y) = x; 


u 2 = 0 * 


(1-16) 
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These two modes represent the motion of the x-y plane while the re 
maining modes signify elastic motion with respect to the x-y plane. 
Hence these elastic deformation modes are identical to those in the 
case of free vibration of the aircraft, and satisfy Equations (1-1) , 
(1-2) and (1-4). With the family of normalized mode shapes, so speci- 
fied, the equations of motion can be simplified to 


ME + id M E = Z 
n s n n n^n n 


n = 1,2. 


(1-17) 


where 

M n = fj*. n pdxdY 

fch 

is the generalized mass of the n mode and 


(1-18) 


Z n = ff s Vn dxd * (1 " 19) 

'til 

is the generalized force of the n mode. 

As previously discussed, the particular values of and depend 

upon the normalization scheme employed. For example, if Equation 

(1-14) were used then we would have M n = for all n. The method of 

normalization has no effect whatsoever on the calculated motion, for the 

product <I> N S N is independent of the particular scaling value used to 

fix d> . This is easily seen from Equation (1-17) together with the 
T n 

definitions of generalized mass and the generalized force. 

A further note on the simplifications in the previous analysis 
is due. The free aircraft mode shapes derived from Equation (1-7) 
are slightly different than the manufacturer's mode shape data for a 
specific flight condition, such as those presented in Figure 2. This 
is because the modal excitations will produce small changes in local 
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angle of attack producing velocity dependent aerodynamic damping 
terms. The net effect of these terms is to change slightly the natural 
mode shapes and natural frequencies of the airplane. Similarly, the 
equations governing the normal coordinates (1-17) ignore structural 
damping effects. These terms are generally accounted for by the 
introduction of a small constant damping term (0.01 to 0.03) in the 
differential equations for It should be clear that the introduction 

of these additional effects may couple the modes of the unrestrained 
aircraft. 

B. Rigid Body Motion 

The preceding analysis of the forced motion of the elastic air- 
frame leads directly to equations governing the pitching and plunging 
motion of the aircraft, (A. 5-11) and (A. 5-10) respectively. However, 
these equations are derived without consideration of aerodynamic 
effects, which although negligible for small elastic deformations, may 
be significant for larger rigid body motions. Hence it would seem 
beneficial to derive expressions for the rigid body motion of the 
airplane incorporating the aerodynamic forces acting on the structure. 

It should be clearly understood that a reformulation of the equations 
governing the rigid body motion will not affect the validity of Equation 
(1-17) , governing the elastic deformations. The derivation of Eq. 

(1-17) is in no way dependent upon the particular time functions for 
translation and rotation, ^(t), and (t) , [See Appendix A, Section 
5], but rather upon the mode shapes for these motions, which remain 
unchanged. 

The rigid body motion is most easily examined by a small pertur- 
bation analysis about a level equilibrium flight condition. The 
resulting ordinary differential equations are not cumbersome, and 
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are quite adequate when the deviations from nominal flight conditions 
are small. We are concerned with the pitching and plunging motion of 
the aircraft described by the following set of equations: 


^0 u _ u 

qS U Q x u U 0 


C a + ^ 0 = 0 
x a qS 


( 1 - 20 ) 







C 6 


z 6 


( 1 - 21 ) 



u 








m. 


( 1 - 22 ) 



(1-23) 


where the perturbation quantities are: u, velocity; a, angle of attack; 
6, pitch angle; h, altitude; and 6, elevator deflection. The equations 
(1-20) through (1-23) are examined and put in the form of a minimum 
realization linear system in Appendix B. The interested reader is 
directed to Reference 5 for a complete derivation and a more exhaustive 
study. 

The effect of the rigid body dynamics is to superimpose an oscil- 
latory motion on the constant altitude nominal flight trajectory. This 
motion, the short period mode of the airplane is characterized by a 
low frequency sinusoidal impulse response in angle of attack and pitch 
rate perturbations. The center of gravity location is most important 
in the short period mode response as it determines the aerodynamic 
restoring moments due to a perturbation in angle of attack. The 
further forward the center of gravity, the larger the restoring 
moments . 
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Approximate dynamics for this motion are obtained by examining 
the homogeneous response of the dominant second order system in angle 
of attack and pitch rate perturbations: 




mU, 


qS 


0 


(1-24) 


m 


a + 


qSc 



0 


(1-25) 


Disregarding the term as it is most often negligible, the charac- 

a 

teristic equation for this motion is easily shown to be 


X 2 + 


2dw X 
n 


w 2 = 0 
n 


(1-26) 


where 


and 


..2 _ z a m 6 2U 0 

m 

a 

n mu. 

i 

I 

0 

y 

y 

qS 

qSc 

qSc 


* C 

? — c 


z 

2U_ m: 

2dw — — 

a 

, o 0 

n 

mU. 

+ I 


0 

_y 


qS 

qSc 


(1-27) 


(1-28) 


In the SST model, the center of gravity is quite far aft and the wing 
loading is light, resulting in small restoring moments and a relatively 
long period for the rigid body motion. The approximate values for 
natural frequency, and damping ratio, d, are calculated to be 

1.48 sec ^ and 0.177 respectively. 

C . The Control 

For this investigation, the control surfaces will be restricted 
to the elevators. In the swept wing B-2707 the elevators represent 
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the major control surfaces with which the pilot can vary forces in the 
z-direction. As would be expected, elevator deflections are an efficient 
means of producing rigid body motion of the aircraft, particularly 
pitching moments. Furthermore, as the elevators do not coincide with 
the nodes of any significant elastic deformation modes, these control 
surfaces should be effective in both stimulating and controlling 
aircraft flexure. 

In the sweptwing configuration, the wing and tail form one con- 
tinuous lifting surface. (See Figure 1). If the elevator deflection 
is non-zero, i.e. the elevator and the wing surface intersect at a 
non-zero angle, then the elevator acts as a new lifting surface with 
an angle of attack equal to the elevator deflection, measured downward 
as positive. From elementary supersonic flow theory, the coefficient 
of lift of a flat plate at an angle of attack a is 

C L = 4a(Mj - 1)~ 1/2 (1-29) 

where M Q is the Mach number. Hence the normal force per unit area of 
the elevator, per unit elevator deflection is given by 

F z = 4(M 2 - l) _1/2 q (1-30) 

where q is the dynamic pressure (one half the atmospheric density 
times the square of the aircraft velocity). 

D. The Cost Criterion 

As mentioned in the introduction, the control objective of this 
study is to minimize the passenger discomfort associated with the rigid 
body and flexing motion of the airplane due to wind gusts. This 
discomfort arises from two quarters. First there is the effect of 
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undesirable accelerations of the aircraft structure yielding a rough 
ride. Second is the low frequency roller-coaster type of motion which 
even without significant accompanying acceleration, can provide a most 
unpleasant flight experience. The discomfort is the result of a number 
of factors including the pitch motion of the aircraft, changes in sound 
level of the engines and low level vertical accelerations, all occuring 
at the same low frequency. Although it is difficult to assess all of 
these factors, an overall measure of discomfort is the magnitude of 
changes in altitude that result from the roller coaster type motion. 

Studies of human vibration response have shown that human dis- 
comfort varies not only with the magnitude of disturbing accelerations, 
but with their frequency as well. People seem to have an extreme 
distaste for accelerations applied at frequencies between four and 
eight cycles per second. This particular phenomenon is graphically 
shown in Figure 4, which represents tests conducted in concert with 
SST flexure studies. 

Figure 4 suggests that passengers riding the SST would be fairly 
sensitive to the accelerations induced by the aircraft flexure, occurring 
at frequencies from 1.47 to 3.81 cps. We might then expect if the gross 
motion due to flexure were small enough, that the primary influence 
of the flexure modes on the ride comfort is the acceleration environ- 
ment of the passenger compartment. This is exactly the situation 
that has occurred in flight simulations. Although the elastic motions 
were considerably larger than for subsonic airliners, they were still 
small enough so that the acceleration environment was the predominant 
factor. [Ref. 2, p. 84] From Figure 4, we see that the low frequency 
rigid body accelerations (the dominant rigid body frequency is 0.232 
cps) need to be quite large compared with the flexure accelerations 
before they exhibit the same degree of unpleasantness. However, the 
gross vertical motion of the airframe due to pitching and plunging may 
be orders of magnitude larger than the elastic motion. Thus it is 
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likely that except for conditions of high turbulence, (which in genera] 
will occur at lower altitudes) the vertical motion due to the rigid 
body dynamics rather than the rigid body accelerations will have the 
predominant effect on passenger discomfort. 

A reasonable control objective would thus be to minimize the 
square acceleration caused by the airplane flexure and the square 
perturbation in altitude caused by the rigid body motion over the dura- 
tion of cruise. The critical independence of the response modes of 
the unrestrained airplane allows us to decouple the total motion of 
the aircraft, and weight separately characteristics of the disturbing 
behavior. The orthogonality of mode shapes will lead to considerable 
simplification of the functional representation of the control objective. 

Let us deal first with the accelerations effected by the airplane 
flexure modes. From equation (1-17) we have 

= ““n C n + H- ff r z*n dxdy • (1 ' 31) 

n S 

From Section C we know that the normal force is restricted to the 
elevator and is proportional to the elevator deflection. Hence, in 
light of Eq. (1-30) we can write 

4 = -“n^n + ( W 6 d- 32 ) 

where 

G n = X£ 4(M 0 " 1 )"" 1//2 q<!> n (x,y)dxdy (1-33) 

and the integration is restricted to the surface of the elevator. 
Multiplying Equation (1 — 32) by <j> n and then summing over the deformation 
modes (n = 3,4,...) yields the total acceleration due to flexure 
dynamics : 
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(1-34) 


w p (x,y,t) = l <j> n (x # y) [-“^ n ( t> + FT 6] 


As the acceleration of the fuselage is of greater import than the 
accelerations of the lifting surfaces, we would like to weight the 
fuselage more heavily than the wings and tail when integrating the 
squared acceleration function over the surface of the airplane. While 
arbitrary weighting functions may allow considerable latitude in the 
formulation of the cost functional, one particular function affords 
considerable simplification in the form of cost. As the mode shapes 
are orthogonal with respect to the mass distribution p(x,y) [See 
Appendix A, Section 4] use of this as a weighting function will mean the 
absence of cross products of modes in the cost functional. In addi- 
tion, for the idealized SST model, the mass distribution is signifi- 
cantly higher over the fuselage than for the lifting surfaces, and 
thus readily satisfies the original condition of a priority weighting. 
Accelerations of large mass elements are thus heavily penalized by 
this weighting scheme. It is reasonable to expect that the resulting 
control system will not allow large accelerations of massive structural 
subassemblies and hence will not produce excessive structural loading 
to achieve ride quality. 

From Equation (1-34), we have the "cost" of the disturbing flexure 
induced accelerations, where T is the length of the cruise, as 


PdXd ydt 


r t r r 

, G 

fff 

l *„<-»» V H- S » 

o s 

_n n 


pdxdydt. (1-35) 


Let us denote the left hand side of Equation (1-35) by J p . Expanding 
the bracketed term in the equation and making use of the orthogonality 
of the mode shapes we can write 
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J F = / ^7X^n pdxdy[ “n 5 n 


2<*l ^ E 6 + (^-) 2 6 2 ]dt 

n M n M 

n n 


(1-36 ) 


which simplifies to 


/•T . ~ ^ G G 0 ~ 

J = / y M [co 4 C 2 - 2 oj 2 £ « + (^r— ) 2 6 2 ] dt 

F J n L n n^n n M M 


(1-37) 


For the rigid body inodes, although we want to minimize vertical 
motion rather than accelerations, the approach is quite similar. The 
gross rigid body motion is given by 


w R (x,y,t) = h + x0 • 


(1-38) 


Squaring w_, weighting by the mass distribution and integrating 

X\ 

over the surface of the airplane and the duration of the flight yields 


/ 0 j(jpR pdxdydt = 


+ 2xh6 + x 2 0 2 ) p (x,y)dxdydt .(1-39) 


Letting J D denote the left hand side of (1-39) , we can expand the 
integral on the right to obtain simply 


/' T 2 /" T 2 

J n = M. / h dt + M _ / 0 dt • 

R IJq ZJ Q 

The entire cost for input disturbances then becomes 


(1-40) 


J = C F J F + C R J R (1 " 41) 

where the weighting coefficients c„ and c„ can be adjusted to obtain 

r K 

any balance of control effort between flexing and rigid body motions. 
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The reader should notice that although the cost functional for 
the fuselage accelerations. Equation (1-37), does not include the 
damping terms mentioned in Section A. 3, there is no difficulty in 
extending the cost to include this feature. Thus, if we add to the 
right side of Eq. (1-32) a term representing damping, -2< 3 n “ n ^ n ' where 
d^ stands for the damping ratio of the nth mode, Jp becomes 


r T 

= / I 

•'A TO 


M [tu 4 C 2 + 4d 1 + 

n n n n n n n 


aa 2 2 i2 
4d a) £ 
n n n 


2 “* ir 5„ 8 

n 


G G 9 9 

- 4d a £ 6 + (w^) 0 Jdt 
n n M M 

n n 


(1-42) 
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CHAPTER 2 

T he Control Problem 

This chapter of the thesis deals with the formulation, solution 
and implementation of the mathematical control problem motivated by 
the analysis of Chapter One. Section A presents the framework and 
solution of the problem. The computer implementation is considered 
in Section B. As in the preceeding chapter, detailed calculations have 
been relegated to the appendices. 

A. Solution of the Control Problem 

From the modeling of the physical system in Chapter 1, it is 
clear that the dynamics of the rigid body motion and the first four 
elastic modes can be combined in the framework of linear ordinary 
differential equations. Let us define the thirteen (13) component 
state vector x as follows: 

x' = [u a e 0 h K 3 c 4 ? 4 ?5 ?5 S 6 ^6 1 * (2_1) 

(To avoid confusion between vectors, matrices and scalars, the notation 
of lower case underlined letters for vectors, upper case underlined 
letters for matrices will be used. Scalars will not be underlined. 

The prime superscript on vectors or matrices denotes transpose.) 

Note that x (t) is a differentiable function of time as are all its 
components . 

We will define the matrix A as below: 
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( 2 - 2 ) 


where A^, the system matrix of the rigid body motion,, is given by 
[See Appendix B ] 


C x 
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C U./a 
x O' 


C z /a 
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C z /a 
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2U 0 V / *’ 


b{-C + ^ C C /a) 
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0 b == (C .+ C ) 
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0 a 0 


(2-3) 
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with 


raU 0 

a = , and 

qS 


b = 2|£ 

y 


The matrices of the elastic motion are defined 
~ 0 1 


A. = 
-3 


2 

-a) . 
• 3 


-2d .(i) . 
3 3 


j = 3, 4, 5, 6 


and the 0 signify matrices all of whose values are zero. 

Continuing, we define the control sensitivity matrix B as 

B. . ,0 C.A 0 C^/a, 0 0 

G 3 /M 3 0 G 4 /M 4 0 V M 5 0 W ' 

It is easily seen that Equations (1-17) , (1-20) , (1-21) , (1-22) , 

(1-23) can be combined in the matrix form 

x ( t ) = Ax (t) + B6 (t) • 

A schematic representation of the system is shown in Figure 5. 



x(t) 


Figure 5: Block Diagram of System of Equation (2-6). 


(2-4) 


(2-5) 

and 

( 2 - 6 ) 
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The cost functional can similarly be put into matrix notation in- 
volving the integral of a quadratic form. Let us construct the 
13 x 13 matrix L as follows: 

L 3 , 3 = c R M 2 

L 5 , 5 = C R M 1 (2-7) 

L 2j,2j = C F M j“j j = 3 ' 4 ' 5 ' 6 

and all other components of L are zero. We define the 12 x 1 matrix 
N as 


N' = [0 


-C F W 3 G 3 


" c F u 4 G 4 


- c F u 4 G 5 


■ c F w 4 G 6 


0 ] 


( 2 - 8 ) 


Finally construct the lxl matrix R 


6 G 

5 = l vr 

n=3 n 


(2-9) 


With the matrices L, N, and R defined as above, the cost functional 
given in Equation (1-41) becomes 



( 2 - 10 ) 


and the control problem can be fully stated as follows: minimize J 
subject to Eq. (2-6) . 
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A solution for the problem exists provided that the matrix L is non- 
negative definite, R is positive definite, and L - NR 1 N' is non- 
negative definite. Equations (2-7), (2-8), and (2-9) substantiate 

that these conditions are met. The minimizing control 6*(t) is given 
in the form of a feedback control law. 

6 * (t) = -R _1 (N'+B'K(t) )x(t) (2-11) 

where K(t) satisfies the matrix Riccati equation 

K (t) = -K(t)A - A'K(t) + (K(t)B + N)R _1 (N' + B'K(t))-L; 

K (T ) = 0 . (2-12) 

Furthermore, the minimum cost is given by 

x' (0)K(0)x(0) . (2-13) 

The reader is directed to Appendix D for the derivation of the optimal 
control . 

A few words on the meaning of the control law are in order. The 
control law. Equation (2-11) is independent of initial values of the 
state. Hence 6* is the minimum cost control for arbitrary initial 
conditions on x(t) which have been for this reason left unspecified 
in Eq. (2-6). As T, representing the time of cruise is large compared 
to response times of the system, the effect of the control is obviously 
to drive the state towards zero. With 6* specified by Eq. (2-11), we can 
substitute into Eq. (2-6) to achieve the minimum cost trajectory 

x*(t) = (A - BR -1 (N 1 + B'K(t) )x* (t) . (2-14) 
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Note that the control 6* is memoryless, hence at any time t in the 
interval from 0 to T , we can consider the problem as minimizing the 
cost from t to T with initial state x(t). 

So far the equations of motion for the airplane have neglected the 
input disturbances, random wind gusts. Let us now consider this effect 
and rewrite Eq. (2-6) as follows: 


x(t) = Ax (t) + B6(t) + Dg (t) 


(2-15) 


where g(t) serves to model the random wind gusts, and D is the in- 
f luence matrix for these gusts. The wind disturbances g(t) are 
clearly independent of the state, x(t) , and the control & (t) . Further- 
more, we may assume that the correlation times of g(t) are short 
compared with characteristic times of the system. This last assump- 
tion implies that knowledge of the gust history up to any time t, 
does not allow for effective prediction of future gust inputs far 
enough in advance to influence the control policy. 

In terms of the stochastic differential equation of state Eq. (2-15) , 
it is no longer valid to minimize the deterministic cost J. Rather, 
the control objective must now be to minimize the average or expected 
value of the cost 



( 2 - 16 ) 


An important result of stochastic control theory states that if 
the input disturbances are white noise, as assumed above, then the 
control to minimize J, given by Eq . ( 2 -ll) is identical to the mini- 
mizing control for J subject to the stochastic system of Eq. (2-15) 
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The reader is directed to Reference 4, Chapter 14 for a more complete 
development. 

The system. Equation (2-6) or (2-15), is time-invariant as the 
matrices A, B, and D, are constant, (the disturbance process, g(t) 
is assumed stationary) and all time varying signals will become very 
nearly ergodic shortly after the initial time. As the length of cruise, 
T, is extremely large compared to characteristic response times for the 
airplane, we would expect little change in either the total cost or 
the minimum cost trajectory for an increase in T. As T grows larger, 
the feedback gain matrix K(t) tends to a constant. The minimizing 
control is still given by Equation (2-11) , but the matrix K is 
now given by the algebraic equation 

-KA - A'K + (KB + N)R -1 (N' + B'K) - L = 0. (2-17) 

Sufficient conditions for a well posed solution, i.e. a positive 
definite solution to the algebraic equation (2-17) and a finite 
minimum cost, include that the system 



Figure 6. Block Diagram of System of Equation (2-15) with 
Closed Loop Minimum Cost Control 
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x = Ax + Bu ; y = Cx 


(2-18) 


where C is the square root matrix of L, i.e. 


L = CC' (2-19) 

(the existence of C is guaranteed for L symmetric and non-negative 
definite) be completely controllable and completely observable. The 
uncontrolled system need not be stable, yet the feedback controlled 
trajectory will be asymptotically stable. The SST model satisfies 
these conditions. The reader is directed to Reference 3 for further 
details. 

B. Computer Implementation 

The use of computers was restricted to two distinct areas for 
this study, data processing and simulation. The minimum cost control. 
Equation (2-11) , was calculated by digital computer for different 
flight configurations. These results were then incorporated into 
an analog computer simulation where the effectiveness of the control 
could be observed. 

The main thrust of the digital computer solution of the control 
problem lies in the solution of the time-invariant algebraic Riccati 
equation (2-17) . The standard method is a backwards in time integration 
of the differential equation (2-12) until a steady state is achieved. 

For this analysis, a different solution technique was used based upon 
an algebraic algorithm developed by Potter. [Ref. 8] . The method 
devised by Potter involves calculating the eigenvalues and eigenvectors 
of the matrix V, 
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V 


( 2 - 20 ) 


(A - BR _1 N' L - NR~ 1 N' 

BR~ 1 B' - (A - BR -1 N' 

The solution of Equation (2-17) , which may be easily rewritten as 

0 = -K (A - BR _1 N' ) - (A - BR -1 N')'K + KBR~ 1 B’K - (L - NR~ 1 N 1 ) 

( 2 - 21 ) 

can be formulated in terms of the eigenvectors of V. Furthermore, the 
eigenvalues of V can be shown to be the eigenvalues of the matrix 
A - br' 1 !!) 1 + B'K) i.e., the poles of the minimum cost closed loop 
system, Equation (2-14). These assertions are proved in Appendix 
E. 

Potter's method, being an algebraic rather than an iterative 
scheme, has several advantages over the integration methods. The 
user need not be bothered choosing an integration step time, which 
may be crucial for the convergence of the solution. Also, the method 
is less sensitive to initial errors which tend to propagate in an 
iterative scheme. However a problem common to all digital calculations, 
limited precision, can have a devastating effect on the solution 
accuracy, and care must be taken to minimize the range of magnitudes 
of the input quantities through appropriate scaling. 

Other subprograms necessary for the digital calculations included 
MAIN, an input-output routine to set up the matrices in the Riccati 
equation, and then calculate the feedback control gains; POTTER, 
which solves the Riccati equation; an eigenvalue-eigenvector routine 
called ALLMAT (authored by J. Rinzel and R. Funderlic of Union Carbide 
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Corp; Oak Ridge, Term.); and a matrix inversion subprogram TDINVR. 

The analog simulation used the results of the digital calculations 
to examine the behavior of the controlled and uncontrolled systems. 

Due to the limited range of the analog computer, a number of extremely 
small parameters were set to zero. This, however, represented little 
loss of generality as the deleted quantities were orders of magnitude 
smaller than accompanying signals, and their removal ameliorated 
the problem of signals being swamped by noise. In addition to moni- 
toring the effects of the control, the analog computer was a useful 
tool in assuring that the data inputs for the model were representative 
of predicted SST performance criteria. 
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CHAPTER 3 


Application to the Model and Results 

In this chapter the theories explained in the previous chapters 
will be applied to an idealized model of the B-2707 prototype. The 
first segment will deal with the physical aspects of the test airplane. 
Results of the digital computations and analog simulation will be pre- 
sented in later sections. Remarks are included in the final section. 

A. Application to the Model 

A. 1. Generation of the Test Model 

In preparing the analysis for computer implementation, it became 
necessary to create an idealized version (Fig. 1) of the B-2707 pro- 
totype aircraft. This model differs very slightly from the actual aircraft 
in overall size, but as data on the mass distribution of the B-2707 
was not available, it was necessary to postulate a distribution func- 
tion for the model. (It should be clear that the problem formulation 
is not dependent upon the specific mass distribution of the aircraft; 
the latter enters solely as a data input and even then, smoothed 
via integration.) A two phase mass distribution was chosen (constant 
over the fuselage surface and constant over the lifting surfaces) as 
the simplest form which would allow specification of two critical 
parameters, the aircraft mass and center of gravity. These were chosen 
to coincide with appropriate values for the actual B-2707, 11,700 slugs 
and body station (B.S.) 2400 inches respectively. 

To facilitate calculation of aircraft properties, the model was 
partitioned by the superposition of a rectilinear grid into approximately 
seventy-five segments. The extreme forward and aft fuselage sections. 
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comprising the first five hundred (500) inches and last two hundred 
(200) inches of the body were excluded from further consideration. 

These fuselage extensions are uninhabited (the crew compartment begins 
at B.S. 560) and vibration mode data was unavailable for them. Modal 
deflection was determined by averaging the deflection at the four lattice 
points of each segment, and was assumed constant over the segment. 

The generalized masses M n and generalized force coefficients, 

C>n (generalized force per unit elevator deflection) were calculated for 
the flexure modes by integrating over the grid of the aircraft model. 

The manufacturer's mode shapes shown in Figure 2 were used as the 
normalized deflection modes for these integrations. The two sets of 
parameters, and G n , together with the modal frequencies were pre- 
calculated and fed into the digital programs as input. In addition, 
to compensate for the structural damping present in the flexure dynamics, 
but ignored in the derivation of Equation (1-17) , a term corresponding 
to a damping ratio of 0.03 was included in each of the equations of 
flexure motion. 

A. 2 Flight Environment and Parameters 

For the purposes of this investigation, the aircraft was assumed 

to be in steady state cruise. The altitude and mach number of cruise 

were chosen to be representative of proposed SST operating conditions, 

65, -000 feet, and Mach 2.7 respectively. As the speed of sound at this 

altitude is 967 feet per second, the aircraft velocity is 2600 feet 

per second. The mean air density at 65kft is 1.78 x 10 -4 slugs per 

cubic foot yielding a dynamic pressure of 602 slugs per foot-second 
2 

squared (lbs/ft ). The reader desiring the calculations of these and 
further quantities is directed to Appendix B. 

The mass distributions over the wings and fuselage were calculated 
to be 0.0042 and 0.0163 slugs/sq. in. respectively, and the lifting 
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surface of the model is computed to be 6680 square feet. The reference 
wing chord was chosen to be 50 feet yielding a pitching moment of 
inertia of 1.875 x 10 7 slug-ft 2 . With these aircraft parameters at 
hand, the lift and drag coefficients are readily calculated as 0.093 
and 0.541 x 10~ 2 . From these values, the equilibrium angle of attack 
becomes 0.058 radians, or approximately 3.3 degrees. 

The stability derivatives appearing as coefficients in Equations 
(1-20) thru (1-23) and examined in Appendix B can be calculated (or 
in some cases approximated) from the parameters determining the flight 
environment above. For the idealized B— 2707 model, these coefficients 
have the values (from the defining equations of Appendix B) : 
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-0.0108 
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= negligible 
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= 
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m 


(3-1) 


Moving for a moment to the flexure motion, the generalized 
masses. Equation (1-18), were computed via integration over the par- 
titioned model. As by convention, we have designated the aircraft 
mass and pitching moment of inertia as and respectively. The 
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remaining generalized masses associated with the flexure modes (mea- 
sured in slugs) are 

[M 3 , M 4 , M 5> Mg] = [1300, 63.9, 411, 726] • (3-2) 

Similarly, the generalized force coefficients. Equation (1-33) , are 
computed via integration over the elevators of the partitioned model 
(See Figure 1), and for this flight condition take on the values, 
(measured in lbs/deg elevator deflection) 

[G 3 , G 4 , G 5 , G g ] = [5290, -1500, -2330, -616]. (3-3) 

It should be noted that the Equations of rigid body dynamics 
(1-20) thru (1-23) and hence the matrix A^ defined by Equation (2-3) , 
are based upon the angle perturbations being measured in radians. 

For computational reasons dealing with scaling, it is advantageous to 
measure angles in degrees. Therefore several entries in the 
matrix need to be adjusted by the radian to degree scale factors, 

57.296. Henceforth, we will use the units common to the computer 
implementation of the model; velocity is feet per second, angles of at- 
tack, pitch and elevator deflection are measured in degrees, altitude is 
in feet and flexure deformations are measured in feet, and flexure fre- 
quencies are radians per second. Designating insignificant or negli- 
gible terms by *, the matrix A^ becomes 

-1 . 43xl0 -3 -0.546 -0.562 0 0 

* - 0.212 0 1 * 

A 1 = 0 0 0 10 (3-4) 

* -2.14 0 -0.31 * 

0 -45.4 45.4 0 0 
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and the remaining non-zero blocks of the A matrix, the matrices of 
elastic motion are: 



0 

1 


0 

1 

l> 

u> 

II 

-85.2 

-0.554 

*4 = 

-200. 

-0.848 


0 

1 


0 

1 

II 

in 

<1 

-259. 

-0.965 

—6 = 

-568. 

-1.43 


The control sensitivity matrix B, defined by Equation (2-5) takes on 
the value 

B' = [0 0.0191 0 -2.03 0 0 4.06 0 

-23.5 0 -5.66 0 -0.848 ] . (3-6) 

The cost, defined in Equation (1-41) has a variable element, the 

ratio of C to C„, which allows for a range of minimum cost controllers 
R t 

to be specified. For the digital computation of the feedback control 
gains, these two parameters were normalized so that C F was identical 
to unity and the weighting coefficient for the costs attributable 
to rigid body motion, denoted as the cost ratio, was defined 

CR = C R /C F . (3_7) 

In this completely equivalent formulation of the quadratic cost func- 
tional, the matrix N, defined in Equation (2-8) becomes 
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M' = [0 0 0 0 0 -4.5 x 10 5 0 3.0 x 10 5 

0 6.02 x 10 5 0 3.50 x 10 5 0]. (3-8) 

The diagonal of the L matrix, Eq. (2-7), takes on the values 

L = [0 0 5.71 x 10 3 CR 0 1.17 x 10 4 CR 9.45 x 10 6 

0 2.55 x 10 6 0 2.75 x 10 6 0 2.34 x 10 8 

0 ] • (3-9) 

Finally, the R matrix, defined in Equation (2-10) becomes 

R = 7.05 x 10 4 . (3-10 

For reasons of computational accuracy, the elements of the N, L, 
and R matrices were further scaled down by a factor of 10 8 in the 
final digital programs. 


B. Results 


B. 1 . Digital Calculations 

Four different values of CR were used in the generation of digital 
solutions to the control problem in an attempt to view the effect of 
this parameter on the behavior of the closed loop feedback system. As 


CR varied from 0.05 to 100., with intermediate values of 1., and 10., 
the relative weightings assigned to rigid body and flexure perturbations 
varied by a factor of 2000. The particular values for CR were chosen 
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arbitrarily and in no way represent an exhaustive set or range of cost 
strategies. Furthermore, it should be clear that each value of CR 
leads to its own minimum cost solution. By varying CR we are not 
searching for an extremum again, but rather only viewing the consequences 
of a set of minimum cost trajectories. 

The major computational effort of the digital program lies in 
the solution of the algebraic Riccati equation. Equation (2-17) , 
which is used to generate the minimum cost feedback gains 

fb = -R _1 (N’ + B'K) (3-11) 

presented in Table 1. As CR increases, the rigid body influence on 

the control grows, as evidenced by the increases in fb a , fb 0 , » 

and fb. . However, of the flexure motion, only the feedback gain for 
h 

the first mode shows an appreciable decrease. The gains for the 
third and fourth modes instead increase with CR. This situation 
is counterintuitive at least, and helps to point out some of the 
subtleties encountered in generating minimum cost feedback controllers. 

As stated in Chapter 2 and in Appendix D, the solution of the 
algebraic Riccati equation also yields the poles of the closed loop 
minimum cost trajectory. These pole locations have been plotted 
in Figure 7. The reader is cautioned that this diagram should not 
be interpreted as a root locus. The points plotted represent the 
poles of five distinct systems, not merely one system with a changing 
loop gain. Although only one parameter, CR, is varying in four of 
the systems, this parameter is involved in the closed loop system func- 
tion in a far from trivial way. Secondly, and less profound, the 
reader should note the differences in scale between the Re(s) and Im(s) 
axes, chosen to help illustrate the pole locations. 
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6.72 


2.66 

fb r 

^ 3 

= 

-0.178 

/ 

-0.573 

f 

-0.872 

f 

-1.05 

fb 

^ 4 


-2.64 


-1.64 


-0.607 


0.724 

fb- 


0.301 


0.334 


0.350 


0.349 



3.69 


5.92 


7.81 


9.93 

“t. 


0.780 


0.723 


0.664 


0.515 

fb C 6 


10.7 


13.1 


15.2 


17.5 

K.- 


0.845 


0.781 


0.708 


0.605 


TABLE 1. Feedback Gains 


As the uncontrolled system is block diagonalizable , [See Equation 
(2-2)] we can interpret the poles of this system in terms of its paral- 
lel subsystems. Referring to Figure 7, the four poles with imaginary 
points greater than 9 are precisely the frequencies of the four flexure 
modes. The remaining two uncontrolled pole pairs together with the 
real pole at -1.43 x 10 3 corresponding to u, the forward speed 
perturbation, represent the longitudinal group rigid body motion. 

The pole at the origin is but one of a pair related to the two states, 

9 and h, which are integrals of other states and devoid of feedback. 

The complex pole at (-.26, +1.46) represents the short period mode, 
and as noted in Chapter 1 Section B, is the dominant second order 
system in angle of attack and pitch rate. [See Equations (1-24 ), (1-25 )] . 
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1. A cluster of symbols indicates that 
the true values lie too closely to- 
gether to be discerned at this scale, 
and should be taken coincidend at the 
center of the cluster. 

2. All five systems exhibit a real pole 
at the point (-1 . 43*10~3 , 0) , not 
shown on this graph. 
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B. 2 Analog Output 

The open loop system model of Equation (2-6), was simulated on an 
analog computer. As the computer was of limited voltage range, con- 
siderable scaling was necessary before the simulation could be imple- 
mented. The negligible terms of the block system matrix A , Equation 
(3-4) , were set to zero as their slight magnitude precluded precise 
representation on the machine. The block diagram for the simulated 
system is shown in Figure 8. The feedback gains from the solution of 
the control problem were incorporated into the simulation for the 
CR = 1 and CR = 100 trials, yielding models of the minimum cost feed- 
back controlled trajectory. Equation (2-14), for these cases. 

Flight disturbances were simulated in a number of ways. First 
a standard elevator pulse of magnitude 5° for a duration of one second 
was used as a test input. This input is effective as a means of ex- 
citing both rigid and flexible body modes and demonstrates the ability 
of the autopilot to damp these modes. However, the elevator pulse is 
not representative of gust inputs and is not effective in demonstrating 
the trade-off of flexible body acceleration and rigid body displacement. 
Hence, responses were obtained for a standard set of initial conditions 
representing an initial nose up rotation of the aircraft by 5°. The 
resulting responses for CR = 1 and CR = 100 vividly demonstrate design 
trade-offs. Finally a vertical gust represented by the waveform 
1-cos wt was used to illustrate the design trade-off as it would apply 
to the gust input case. 

We consider the elevator pulse input first. It is represented as 
6 0 (t) = 5 (u_ 1 (t) - u_ 1 (t-l)) (3-12) 

where u_^(t“tg) is the unit step at t^. As the open loop system 
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is composed of uncoupled subsystems, the responses to com- 
pletely illustrate the pole locations of Figure 7. The an- 
gle of attack perturbation represents the short period mode, 
a damped sinusoid of approximately 1.50 radians per second, decaying 
to zero, with a peak magnitude of five degrees. The pitch angle per- 
turbation is similar in form and magnitude, but as the integral of the 
short period mode, it does not reach a final value of zero. Hence the 
net effect of the disturbance is to rotate the aircraft in space. As 
the glidepath angle (0-a) has a negative final value, the altitude 
decreases steadily. The 0 and a time histories are shown in Figures 
9 and 10. (For all analog computer output to follow the time scale will 
be one horizontal division per second; the vertical scale will be 
individually indicated. In all cases, the input disturbance to the 
system has been 6-(t).) 

The responses of the normal coordinates of the flexure motion 
(when the meaning is clear we will not make any distinction between 
the normal coordinates and the flexure motion they represent) are 
shown in Figures 11 thru 14. As the equation of motion of each 
flexure mode (Equation (1-17)) represents an uncoupled oscillator, 
the flexure responses are damped sinusoids corresponding to the un- 
controlled system poles of Figure 7. The magnitude of each response 
is directly proportional to the corresponding control gain, G n / M n 
[See Equations (2-5) and (3-6)]. The large acceleration for the 
second mode is restricted to the extremes of the lifting surfaces 
[refer to Figure 2] , and is due to the low generalized mass of this 
mode . 

The pole locations of the minimum cost feedback controlled 
trajectories graphed in Figure 7 indicate roughly the system responses 
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to be expected from the analog simulations. In both the CR = 1 and 
CR = 100 trials the poles corresponding to the three high frequency 
flexing modes are coincident, hence we would expect the responses of 
C 4 / £ 5 / and £g to be nearly similar from one simulation to the next. 

For the first mode however, the change in pole location with increasing 
CR indicates that the fundamental frequency of this response will also 
increase, with the damping ratio remaining approximately constant. 
Similarly, the frequency content of the short period and altitude 
responses increases with CR. As CR gets larger, from Table 1 we see 
that the feedback gains grow, hence the total loop gain increases and 
responses tend to be faster. 

The pitch angle history, 0 (t) is shown for the CR = 1 simulation 
in Figure 15. The pitch angle motion frequency along with that of 
the angle of attack, Figure 16, is approximately three radians per 
second, and considerably speeded up from the open loop short period 
mode. The responses for 0 and a are nearly identical as h(t), pro- 
portional to the integral of their difference, is being quadratically 
penalized. The rotations reach a peak magnitude just under four 
degrees while the altitude perturbation. Figure 17, remains less than 
twelve feet at its greatest. 

The control history is shown in Figure 18. (Actually - 6 (t) is 
presented in the figures; when -<5{t) is positive, the elevator deflec- 
tion is upwards). The discontinuities of 6 Q ( t ) are evident on the 
control history graph as the spikes at the initial time and one second 
thereafter. The control is primarily sinusoidal in shape with a 
basic frequency of about three radians per second, in close accord 
with the pitch and angle of attack motion. The range of elevator 
deflection for this minimum cost trajectory is large, nearly twelve 
degrees on either side of the equilibrium zero point. (The large. 
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Figure 15. 0 (t) (CR = 1) 
Scale; 0.2 degree/line 
Full scale: 5 degrees 



Figure 16. ct(t) (CR = 1) 
Scale: 0.2 degree/line 
Full scale: 5 degrees 
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Figure 17. h(t) (CR = 1) 
Scale: 1 ft. /line 
Full scale: 25 ft. 



Figure 18. -6 (t) (CR = 1) 
Scale 0.5 degree/line 
Full scale: 12 1/2 degrees 
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control effort is at least partly attributable to the specific form of 
the input.) In addition, a comparison of Figure 18 with Figures 15 and 
16, show that after the external disturbance has been removed, the con- 
trol 6 follows the 0 - a motion, in that maxima of one align with maxi- 
ma of the other. This behavior may be explained in terms of the feed- 
back gains fb 0 and fb a . From Table 1, regardless of CR, these gains 
are always in a ratio of approximately -1.1 to 1. Since the perturbation 
altitude rate is proportional to the flight path angle, or 0 - a, 
the feedback to the controller from these two states can be interpreted 
as feedback of the flight path angle and additional feedback on the 
pitch. As 0 - a is very small, this feedback is roughly 0.1 fb 0 6 and 
for the CR = 1 simulation, corresponds to 2.5 6 . 

The normal coordinates for the flexure accelerations , £3 thru 
j? are presented in Figures 19 thru 22 respectively.* The CR = 1 
controller significantly moves the pole corresponding to the first 
flexure mode, as shown in Figure 7 and documented in Figure 19. The 
peak acceleration for this mode has been reduced by about half, and 
the motion is quite irregular rather than sinusoidal. 

As the second and third modes are quite close in frequency, it 
would be natural to expect that their responses in the closed loop 
system would be similar. This is the case, as evidenced by Figures 
20 and 21. The magnitude of each mode has been reduced by about half 
over the open loop simulation and considerable damping has been added. 
Furthermore, each response has the appearance of the sum of two sinusoids 
a fundamental of frequency about three radians per second and a 
higher frequency component, near fifteen radians per second. Finally, 
the acceleration of the fourth flexure mode indicates little effect 

*(As in the case of the control history, is actually graphed. These 

normal coordinates are abstract quantites and their sign depends solely 
on the orientation of the mode shapes from Figure 2.) 
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Figure 19. -£ (t) (CR = 1) 

J 2 

Scale: 1 ft. /sec. per line 
Full scale: 25 ft. /sec? 



Figure 21. -? 5 (t) (CR = 1) 
Scale: 1 ft. /sec? per line 
Full scale: 25 ft. /sec? 



Figure 20. -i At) (CR = 1) 

’ 2 

Scale: 5 ft. /sec. per lines 
Full scale: 125 ft. /sec? 



Figure 22. g (t) (CR = 1) 
Scale: 0.2 ft. /sec? per line 
Full scale: 5 ft. /sec? 










Figure 27. -£,(t) (CR = 100) 
J 2 

Scale: 1 ft. /sec. per line 
Full scale: 25 ft. /sec? 



Figure 29. 5 ( t) (CR = 100) 

Scale: 1 ft. /sec? per line 
Full scale: 25 ft. /sec? 



Figure 28. -? 4 (t) (CR = 10 
Scale: 5 ft. /sec? per line 
Full scale: 125 ft. /sec? 



Figure 30. g (t) (CR = 100) 
Scale: 0.2 ft. /sec? per line 
Full scale: 5 ft. /sec? 


of the control, except that the magnitude of the response has been 
diminished. 

The system responses for the CR=100 trial agree with what was 
roughly intuited from the pole location map of Figure 7. The pitch 
and angle of attack motion, shown in Figures 23 and 24 is of similar 
form to that from the first simulation, although the fundamental fre- 
quency has increased to roughly four radians per second. More startling 
however, is the reduction in magnitude of the response, which is now 
under 0.6 degrees, nearly an 85% reduction over the CR = 1 simulation. 
Similarly the plunge motion, h(t), shown in Figure 25 is dramatically 
reduced, in this case by an even larger margin. Here too, the basic 
frequency has been slightly increased. The control, 6(t), Figure 26, 
is still aligned with the 8 — a motion, and is of similar, form to the 
control of the CR = 1 trial. Figure 18. However for CR = 100, the 
range of elevator deflection has been significantly diminished, and is 
now limited to five degrees. 

From the plot of closed loop poles from Figure 7, we would not 
expect large changes in the responses of the higher frequency flexure 
modes as CR increases to 100. The graphs of 1^, l 5 , and Figures 

28, 29, and 30 show that this is indeed the case. The responses to the 
second and third modes are not quite as similar as in the CR = 1 trial. 
This is partially explained by the feedback gains fb^ and fb^, which 
for CR = 100 do not act to cancel out the influence of these modes in 
the control. In each of these three responses, the magnitude is close 
to that for the CR = 1 simulation. The first flexure mode response, 
Figure 27, is likewise roughly the same magnitude as for the first 
trial. However, the CR = 100 response is seen to be more regular, and 
of a higher fundamental frequency. For a different input disturbance 
we might expect the first mode to exhibit a larger response for the 


65 





higher cost ratio trial. 

Finally, the sum of the flexure accelerations and the resultant 

displacement of the pitch and plunge motions were examined for both 

simulations at three points along the fuselage, B.S. -800, 

B.S. -2000, and B.S. -3200, corresponding to the forward, middle, and 

rear of the passenger compartment. The acceleration component of the 

n th node at the point (x,y)is given by 0 n (x,y) C n (t). From Figure 

2, we can read off the appropriate value of the coefficient <(i n (x,y) 

and the sum over the modes to achieve the total flexure acceleration. 

The rigid body displacement is computed by adding to the altitude 

perturbation the pitch angle multiplied by the pitching moment arm, 

(x _ -x) . 
c . g • 

The flexure accelerations for the CR = 1 simulation are given 
in Figures 31, 32, and 33. From Figure 2, we would expect the accelera- 
tion at the forward and middle fuselage areas to be dominated by the 
first flexure mode. Comparison of Figures 31 and 32 with Figure 19 
tends to bear this out. The flexure acceleration at the rear of the 
passenger compartment is heavily dependent upon the higher frequency 
modes. At all three points, the accelerations are reduced to about 
half those present in the uncontrolled system, shown in Figures 43, 

44, and 45. 

The CR = 100 simulation yielded similar results for the total 
flexure accelerations, shown in Figures 37, 38, and 39. The magni- 
tudes of the responses were indeed close to those for the lower cost 
ratio trial. Figure 27 illustrates the heavy dependence of the forward 
and middle area acceleration histories on the first flexure mode. 

The rigid body displacements for both trials, CR = 1, presented 
in Figures 34, 35, and 36, and CR = 100, Figures 40, 41 and 42, are 
exceedingly similar in form. The B.S. - 800 response indicates the 
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Figure 31. Flexure acc. at 

B.S. 800 (CR = 1) 

Scale: 0.5 ft/sec? per line 

2 

Full scale: 12.5 ft. /sec. 



Figure 32. Flexure acc. at 

B.S. 2000 (CR = 1 
2 

Scale: 0.5 ft. /sec. per line 

2 

Full scale: 12.5 ft. /sec. 


Figure 33. 

Flexure acc. at 
B.S. 3200 (CR = 1) 

Figure 

34. 

R.B. disp 
B.S. 800 

. at 
(CR • 

Scale: 0.5 

2 

ft. /sec. per line 

Scale : 

0.5 

ft. /line 


Full scale 

: 12.5 ft. /sec? 

Full scale 

: 12.5 ft. 
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Figure 36. R.B. disp. at 
.) B.S. 3200 (CR = 1) 

Scale: 0.5 ft. /line 
Full scale: 12.5 ft. 



Figure 38. Flexure acc. at 
°) B.S. 2000 (CR = 100) 

Scale: 0.5 ft. /sec? per line 
Full scale: 12.5 ft. /sec? 





Figure 39. Flexure acc. at Figure 40. R.B. disp at 

B.S. 3200 (CR = 100) B.S. 800 (CR 

Scale: 0.5 ft. /sec? per line Scale: 0.1 ft. /line 

Full scale: 12.5 ft. /sec? Full scale: 2.5 ft. 



Figure 41. R.B. disp. at 

B.S. 2000 (CR = 100) 
Scale: 0.05 ft. /line 
Full scale: 1.25 ft. 


Figure 42. R.B. disp. at 
B.S. 3200 (CR 
Scale: 0.1 ft. /line 
Full scale: 2.5 ft. 



Figure 43. Flexure acc. at B.S. 800 (uncontrolled) 
Scale: 1 ft./sec.^ per line 
Full scale: 25 ft. /sec? 




Figure 44. Flexure acc. at Figure 45. Flexure acc. at 

B.S. ^2000 (uncontrolled) (B.S. 3200 (uncontrolled) 

Scale: 1 ft. /sec. per line Scale: 1 ft. /sec? per line 

Full scale: 25 ft. /sec? Full scale: 25 ft. /sec? 



effect of the long moment arm (1600 inches) on this motion, while the 
middle fuselage response is almost entirely plunge motion. The B.S. 
-3200 response, with a moment arm of 800 inches is still basically 
the altitude perturbation. However the magnitudes of the response 
are decidedly different. At all three positions, the displacements 
for the CR = 1 simulation are at least three times as large as the 
corresponding displacement for the higher cost ratio trial. 

As previously mentioned, the external input S^ft) is particularly 
successful in exciting the different modes of the aircraft, but less 
effective in contrasting the trade-offs in system performance for 
differing values of cost ratio. The results show a decided decrease 
in the rigid body perturbations for the CR = 100 trial; yet the 
flexure accelerations for the CR = 1 simulations are not appreciably 
smaller. The changes in response for changes in CR are more apparent 
when other means of simulating the input disturbance are used. 

Figures 46 thru 49 illustrate the effect of non-zero initial 
conditions as input to the system. For these outputs, respresenting 
the rigid body displacements and flexure accelerations at the middle 
fuselage area (B.S. 2000), the initial angle of attack perturbation 
was set to five degrees. The initial pitch angle was adjusted to 
approximately 4 . 5 degrees so that the combined feedback from 0 and 
a exactly cancelled out, preventing the immediate saturation of the 
analog computer. The non-zero initial displacement shown in Figures 

46 and 48 is due to the effect of the positive pitch angle. Comparison 
of the results for the two cost ratios more readily indicates the 
trade-offs in system performance. The higher weighting on the rigid 
body excursions results in a reduction by over one half for this motion, 

[See Figures 46 and 48] but the price for this improved dynamics is an 
order of magnitude increase in flexure accelerations . -{See Figures 

47 and 49 ] 
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Figure 46. R.B. disp. at 

B.S. 2000 (CR = 1) 
(from initial conditions) 


Scale: 0.2 feet per line 
Full scale: 5 feet 



Figure 48. R.B. disp. at 

B.S. 2000 (CR = 100) 
(from initial conditions) 
Scale: 0.2 feet per line 
Full scale: 5 feet 


Figure 47. Flexure acc. at 

B.S. 2000 (CR = 1) 

(from initial conditions) 
2 

Scale: 0.2 ft. /sec. per line 

2 

Full scale: 5 ft. /sec. 



Figure 49. Flexure acc. at 

B.S. 2000 (CR = 100) 

(from initial conditions) 
2 

Scale: 0.2 ft. /sec. per line 
Full scale: 5 ft. /sec? 




Figure 54. R.B. disp. at Figure 55. R.B. disp. at 

B.S. 2000 (CR = 1) B.S. 2000 (CR 

Scales 0.2 ft. /line Scale: 0.2 ft. /line 

Full scale: 5 ft. Full scale: 5 ft. 



The characteristics of the specific responses are again at least 
partly attributable to the type of input disturbance employed. The 
initial conditions used were not effective in exciting the elastic 
motion of the airframe, and for the CR = 1 trial did not cause large 
control action. For the CR = 100 simulation, however, the higher body 
feedback gains resulted in a larger control effort which in turn 
excited the flexure modes. 

Finally, a set of runs was made to simulate gust inputs. One 
cycle of the function 1 - cos wt was used to represent a downward gust 
The frequency utilized was 1 cycle per second placing the input near 
resonance of the controlled short period mode. Maximum gust velocity 
was 40 ft. /sec. Figures 50 and 51 illustrate the pitch responses for 
the two cases CR = 1 and CR = 100. We note that as for the other input 
cases, the CR = 100 autopilot damps the pitch motion appreciably better 
than the CR = 1 autopilot. The initial positive pitch in both cases is 
the result of the aircraft weathercocking into the downward gust. 
Figures 52 through 55 show flexible body accelerations and rigid body 
displacements for the CR = 1 and CR = 100 cases at body station 2000. 

As expected, the rigid body motion is better controlled with CR = 100, 
however, the flexible body acceleration is larger by an order of magni- 
tude than for CR = 1. Hence a relatively large penalty is paid, in 
terms of flexible body accelerations for better control of the rigid 
body displacements. Responses at other body stations were similar. 

C. Remarks 

The results of the previous section demonstrate the effect of the 
design parameter CR in the controlled system responses. The feedback 
control in both trials not only diminished the magnitude of the undesir 
able responses, but also lowered the basic frequency of the flexure 
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accelerations [compare Figures 30-33, 37-39, 43-45], which according 
to Figure 4 is a desirable end. If the cost ratio were increased fur- 
ther, we could expect better rigid body responses but a decrease in 
the ability to control the elastic deformations and to lower their 
frequency content. Diminishing CR would lead to the converse behavior. 

No particular value of CR is in an absolute sense best; but rather best 
according to the desires of the designer, who can use this parameter 
to influence the controller according to his own subjective preferences. 

The particular entries of the cost for this investigation, flexure 
accelerations, and rigid body displacements, represent but one of many 
possible variations. There is no difficulty in including the displace- 
ments due to the elastic deformations of the structure, or the rigid 
body accelerations to the cost functional. Similarly there is no 
reason why one could not have additional cost ratio terms, representing 
preferential weightings on other responses. 

Along another vein, the present model could be selectively simpli- 
fied without significant degradation of the responses. The perturba- 
tion in forward speed, u, has been shown to contribute little to the 
model, in that it neither enters the cost nor the control to other but 
negligible degree. Hence Equation (1-20) could be omitted from the 
longitudinal group set and u dropped as a pertinent state variable. 

The fourth flexure mode seems little affected by the controller, and 
since it is of small magnitude, one interested in simplifying the 
model would delete the equations corresponding to this motion. 

In solving the equations of elastic deformation, we made use of 
separability to derive equations for the mode shapes <J> n (x,y), and 
normal coordinates £ n (t). The control problem then used these normal 
coordinates as state variables and developed a control as a linear 
combination of the normal coordinates and other states of the system. 
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If we wish to implement this controller, the rigid body perturbations 
are well understood and easily obtained for feedback purposes. The 
flight path angle, pitch angle, pitch rate, and altitude perturbation 
are available from onboard navigational aids, such as an inertial 
measurement unit. But how do we obtain the normal coordinates for 
feedback purposes? 

A sketch of one method which could be implemented follows. 

Since the feedback controlled system, Equation (2-14) is linear, we 
can transform the system matrix A - BR 1 (N'+B'K) into its Jordan canoni- 
cal form. The resulting transformation will give us a new system 

x (t) = A x(t) (3-13) 

where 

A-BR -1 (N'+B'K) = VAW (3-14) 

and where V and W are orthonormal linear transformations involving the 
eigenvectors of A - BR ^ (N'+B'K), and 


This new system of Equation (3-13) has the property that its state 
variables represent uncoupled subsystems corresponding to the eigen- 
values of the system matrix A - BR 1 (N'+B'K). 

These eigenvalues are exactly those poles plotted in Figure 7. 
Since the elastic motion was assumed band limited, a finite number of 
accelerometers strategically placed along the aircraft would give us 
sufficient data to recover the normal coordinates. The output from 
the accelerometers could be filtered to isolate a frequency corres- 
ponding to one complex set of eigenvalues. The relation between x 
and x then tells us the amount of each state of x contained in a given 
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state of x. Since we can measure the states of x (by integrating 
the filtered output from the accelerometers) , we can thus obtain measure- 
ments for the states of x, i.e. the normal coordinates. 
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APPENDIX A 


Derivation of Equations of Flexure Dynamics 

(The analyses of this section follow from Aeroelasticity , by 
Bisplinghof f , Ashley, and Halfman.) 

A. 1 . Differential Equations of Motion of a Slender Beam 

The equations of motion for free (unforced) vibrations of a long 
slender beam where rotary inertia effects and transverse shear defor- 
mations have been neglected is 

iL (ei + m(x) = o (A. 1-1) 

3x 2 V 3x 2 / 3t 2 

where w(x,t) represents the displacement of the point x along the beam 
at time t, EI is the bending stiffness, and m(x) is the mass distribu- 
tion along the beam. Equation (A. 1-1) is easily derived from a consid- 
eration of force and moment balance on an incremental element of the 
beam. 

Henceforth, we will dispense with this rather cumbersome notation 
for derivatives and will adopt the following scheme: 

0 

f'C) will signify -g— f(') 

£ (•) will signify .-g^r f ( ' ) • 

We note that Eq. (A. 1-1) is separable, i.e. will admit solutions 
of the form 

w ( x , t) = W(x) T(t) . (A. 1-2) 

Substitution of the above into Equation (A. 1-1) yields 

-T = (EI • W" ) (A. 1-3) 

T mW 
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As x and t are independent variables, the quotients in Eq. (A. 1-3) must 

o 

also be independent of x and t. Equating these quotients to to , the 
following ordinary differential equations are obtained: 

T + to 2 T = 0 (A. 1-4) 

(El • W" ) " - ma) 2 W = 0. (A. 1-5) 


Solutions to Eq. (A. 1-5) are eigenfunctions and the values of to are 
eigenvalues of the problem. These eigenfunctions are mutually ortho- 
gonal with respect to the mass distribution, m(x) , as shown below. 

Assume W^, W n , to m , <o n are two eigenfunctions and associated eigen- 
values satisfying Eq. (A. 1-5). Then 


(El *W " ) " - rnto„W m = 0 

m mm 


(El • W ") " - mco W = 0. 

n n n 


(A. 1-6) 
(A. 1-7) 


Multiplying Eq. (A. 1-6) by and Eq. (A. 1-7) by W m and then integrating 
along the beam of length SL we obtain 



(El W ") 
m 

(El w ") 
n 


It 


II 


W dx 
n 


W dx 
m 



mdx 



W mdx. 
m 


(A. 1-8) 
(A. 1-9) 


Subtraction yields 




W mdx 
m n 



[ (El W ") " 
m 


W 


n 


(El W ")" W ] dx. 
n m 


(A. 1-10) 

After successive integrations by parts on the right side of Eq. (A. 1-10) 
we obtain 
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(A. 1-11) 


( a) 2 — a) 2 ) / W W mdx 
m n m n 


{W (EI-W ") ' - W m (El • W " ) ' 
n m m n 


-El (W ' W 1 
n m 


W ' W " ) } 
m n 


l 

O' 


If at least one of the following pairs of boundary conditions holds, 
then the right side of Eq. (A. 1-11) vanishes identically and the 
orthogonality condition is satisfied: 


W = 0 and W 1 =0 

W = 0 and El* W" = 0 
W' = 0 and (El *W") * = 0 

EI-W" = 0 and (EI-W")' = 0. 


(A.l-Ha) 

(A.l-llb) 

(A.l-llc) 

(A.l-lld) 


At this point, ready to tackle the solutions of Eqs. (A. 1-4) and 

(A. 1-5) we shall confine ourselves to uniform beams, hence El and m 

2 

are both assumed constant. Letting a = El/m, the aforementioned 
solutions are easily written by inspection 

T = A sinut + B cosiiit (A. 1-12) 

W = C sinh \J -x + D cosh \/^x + E sin + F cosyP^x. (A. 1-13) 

* a cL 2L « 

where the constants A,B,C,D,E,F are determined by the boundary 
conditions . 

For this investigation, we are interested in unrestrained beams, 
hence the natural boundary conditions are that moments and their 
derivatives with respect to position (shear) vanish at the end points. 
Thus 


W"(0) = W'"(0) = W" (£) = W"'(£) = 0. (A. 1-14) 

Applying these boundary conditions to Eq. (A. 1-13) we obtain 


C = E, D = F 


(A. 1-15) 
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(sinhqi - sinq£)C + (coshqJ- - cosqH)D = 0 

(coshqi - cosq£) C + (sinhqi!. - sinqi.)D = 0 (A. 1-16) 

where q = ^ Eq. (A. 1-16) is a set of linear equations in C and D 
which will have a non-trivial solution only if the matrix of 
coefficients is of less than full rank. This requirement yields the 
following equation in ui, from which the eigenvalues or natural frequen- 
cies of vibration can be ascertained: 


cos ■J —i = (cosh -J — l) 

cl ’’cl 

Solutions to Eq. (A. 1-17) are of the form 


(A. 1-17) 


0 

/ 


for n=0 


ID 


n 


(1.51 2 (j) 2 El/m for n=l 
2 

2n+l ,7T. 2 _ T , _ 

— j- (j) El/m for n=2 , 3 , . . . 


(A. 1-18) 


and the eigenfunctions or normal mode shapes become 

( cosq X--coshq l \ 

i TnEq n A-subq\ J ( 3 inh< 3 n x + sinc J n x ) 

+ (coshq x + cosq x) 


(A. 1-19) 


Ls an arbitrary scaling factor. W (x) has been 

n 


faS 

where q n = and G is 

plotted for n=l,2 in Figure 3. 

It will be advantageous for us to normalize the mode shapes. One 


widely used scheme defines 


<f> (x) = B W (x) 

u n n 


(A. 1-20) 
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where the normalizing constant B n satisfies 


B 


2 

n 


M 



(x) mdx 


(A. 1-21) 


and M is the mass of the beam. Also define M j , the generalized mass 

1.U 

of the j mode as 




(x) mdx 


(A. 1-22) 


Note that for the normalization according to Eq. (A. 1-20) 


M. = M for all j . (A. 1-23) 

1 

With these additional mathematical tools in hand, we are ready to 
proceed to the matter of interest — forced motion of the airframe. 
In this situation Eq. (A. 1-1) becomes 


(EIw" (x, t) ) " + mw(x,t) = F z (x,t). 


(A. 1-24) 


As the normal mode shapes and natural frequencies have already been 
calculated, we can postulate the solution as 


w(x,t) = £ $ (x) 5 n (t) 

n=l 


(A. 1-25) 


where the E (t) are called normal coordinates and remain to be 
n 

determined. Substituting Eq. {A. 1-25) into Eq. (A. 1-24) yields 


00 00 

m.V M n + £ (Eli")" € n = F z (x,t). (A. 1-26) 

n=l n=l 
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Multiplying by <J> m / integrating over the length of the beam, and 
using Eq. (A. 1-8) , we arrive at 


V 1 

n^ 



n m 


mdx 




2 f l f 

E w / <P <P mdx = / 
n nJ Q y n*m jf. 


F z (x,t). 


dx 


m 


(A. 1-27) 


where the order of the integral and summation operators has been inter- 
changed. However the <f> n 1 s are mutually orthogonal with respect to 
the mass distribution, hence Eq. (A. 1-27) simplifies to 


C m + 


w 2 S 

mm 


Z/M 
m m 


(A. 1-28) 


where 


Z 

m 


f l 

J F z (x,t) 4> m (x) dx . 
0 


If F z is a point force of magnitude F z (t) 
we can write 


F z ( x,t) = F z (t) 6(x-x q ) 
and Eq. (A. 1-28) further simplifies to 

C n + w n^n = *n (x 0 > ' F z (t ) / M n • 


(A. 1-29) 


concentrated at x =» x Q , then 


(A. 1-30) 


(A. 1-31) 


^•2. Integral Equations of Motion of a Slender Beam 


For a restrained beam vibrating in the absence of external forces 
we can write 


w(x, t) 



C(x,n) 


w(n, t) 


m(n) 


dri. 


(A. 2-1) 
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where c(x,n) is the conventional influence function measuring the 
displacement at x due to a unit force applied at r, with the origin of 
the x-axis constrained to its initial position. Since Eq. (A. 2-1) is 
separable, it accepts a solution of the form 

w(x,t) = W(x) T ( t) . (A. 2-2) 


Substituting Eq. (A. 2-2) into Eq. (A. 2-1), remembering that x and t 
are independent variables, we obtain 


T + u) T = 0 


' / 
J-l 


w (x) =oj / c(x,n) w(n) m(n) dn. 

-St 


(A. 2-3) 
(A. 2-4) 


Note that Eq. (A. 2-3) is identical to Eq. (A. 1-4), and the family of 
solutions to Eq. (A. 2-4) is identical to the solutions of Eq. (A. 1-5) 
provided that the boundary conditions on Eq. (A. 1-5) and C(x,n) are the 
same . 

Now, let us consider an unrestrained beam, as shown in Figure Al. 
The axes are chosen so that the origin coincides with the center of 
gravity of the beam, the x and z axes align with the principal inertial 
axes of the beam. Further, the beam is assumed symmetric with respect 
to the z-axis. 

Since the beam vibrates freely, the sum of inertial forces in the 
z-direction and inertial moments about the y-axis are zero 



w(x, t) 


dm 



w(x,t) xdm 


n 

T(t) J W 
-St 

fl 

T ( t) J W 


(x) m(x) dx = 0 


(x) m(x) xdx = 0. 


(A. 2-5) 
(A. 2-6) 
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Figure A1 . Beam with Coordinate Axes 


If a unit load in the z-direction is applied at a point r along the 
beam, the resultant acceleration at a point x of the rigid beam is 
given by 



y 


(A. 2-7) 


where M is the total mass of the beam and I y is the moment of inertia 
about the y-axis. The inertial force due to an element of mass m(x)dx 
is thence 


df(x) ~ " < H + “ > <** • (A, 2-8) 

y 

The displacement of the flexible beam with respect to a line 
tangent to the beam at the origin, due to a unit force at the point r 
and the distributed inertial forces, can then be written 

H(x,r) = C(x,r) - Ajtx,*) (I + £_) m(n) dn . (A. 2-9) 

-% y 
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(No attempt will be made to derive influence functions for the beam or 
aircraft. These can be calculated depending upon the physical 
properties of the structure involved. [See Ref. 1, Chapter 2] However, 
in this investigation, little insight into our problem will be garnered 
by pursuing the explicit characterization of the influence function.) 

If beam deflection is measured with respect to the x-axis rather 
than a line parallel to it as in Eq. (A. 2-9), then we can write a 
modified influence function 

G (x, r) = H (x, r) + W(0) + * x (A. 2-10) 


where the last two terms are simply the deflection and slope of the 
beam at the center of gravity. G(x,r) in Eq. (A. 2-10) represents a 
displacement of the beam, hence it must satisfy Eq. (A. 2-5) , yielding 



[H (x, r) 


+ W(0) 


3W (0) 
3x 


x] m(x) 


dx 


0 • 


(A. 2-11) 


However as the origin is the center of gravity of the beam, Eq. (A. 2-11) 
yields directly 


W (0) 


1 

M 


H(x,r) 


m ( x ) dx • 


(A. 2-12) 


Similarly, substituting Eq. (A. 2-11) into Eq. (A. 2-6) and carrying out 
the indicated integration gives 


3W ( 0 ) 
3x 



H(x,r) 


m (x) 


xdx • 


(A. 2-13) 


Having obtained explicit formulae for W(0) 


and 


3W (0) 
3x ' 


we can now 


* We will use this notation to designate the derivative evaluated 
at the point (0) . 
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derive an expression for the influence function of an unrestrained 
beam with respect to its principal inertial axes 


/•£ 

G (x, r) = C (x,r) -/ C(x,n) (i + ^p) m(n) dp 



m(n) dn 


rn+xp 

MI 

y 


— >P)m(n)m(p)dpd£. 

I 

y 


(A. 2-14) 


Analogous to the situation for the restrained beam, Eq. (A. 2-4), 
the integral equation for the freely vibrating unrestrained beam 
becomes 


W(x) 



G(x,r) 


W (r) 


m(r) 


dr 


(A. 2-15) 


Further simplification in Eq. (A. 2-15) is possible. From Eqs. (A. 2-5) 
and (A. 2-6) we see that any terms comprising G(x,r) which include only 
zero or first powers of r will not contribute anything to the integral 
of Eq. (A. 2-15). Examination of Eq. (A. 2-13) shows that the second and 
fourth integrals have this property, and as a result, the influence 
function G(x,r) can be simplified as follows: 


G(x,r) - C(x,r) - / C (n ,r) + ~^)m(n)dn. (A. 2-16) 

y 

Equation (A. 21-5) in conjunction with Eq. (A. 2-16) govern the natural 
vibration characteristics of the unrestrained slender beam. As was 
the case with the restrained beam, the family of solutions to Eq. 

(A. 2-15) is identical to that via the differential equation analysis 
provided that G(x,r) satisfies appropriate boundary conditions. 
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Having derived the influence function of the unrestrained beam, 
we can write the unrestrained beam analog of Eq. (2.2-1). 


r A 

w(x,t) = - / G (x , r) w (r , t) m (r ) dr . (A. 2-17) 

The extension to forced motion of the beam requires some care. Denoting 
the applied force per unit length in the z-direction as F z (x,t), one 
is tempted to write, for forced motion of the beam 


w(x, t) 



G (x, r) 


[-w (r , t) m(r) - F z (r,t)] dr. 


(A. 2-18) 


However, the above is not correct. This is because the notion of the 
displacement function in the free and forced motion is different. In 
the former case, the displacement is measured with respect to the 
principal inertial axes of the beam. Equation (A. 2-15) defines only 
elastic deformation modes of the beam. The effect of the integral term 
in Eq. (A. 2-16) is to subtract off the effect of rigid body translations 
and rotation from the total motion of the beam. Thus Eq. (A. 2-15) 
yields only trivial rigid body (u>=0) motion. For forced motion of the 
beam, the rigid body motion cannot be eliminated to leave the deformation 
modes alone. Thus, in this case, the displacement function, w(x,t) 
must represent the total displacememt including translations, rotations, 

- and elastic deformations. 

The mode shapes generated by Eq. (A. 2-15) are orthogonal with 
respect to the mass distribution, i.e. , 



(x) 


W (x) 
m 


m(x) 


dx 


0 


(A. 2-19) 


for m j- n. Proof of the orthogonality condition for the two dimension- 
al case is given in Section 4. The proof of Eq. (A. 2-19) is exactly 
the same and is not repeated here. 
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The analysis of the forced motion relies heavily on the concept 
of orthogonal mode shapes to represent the motion of the beam, and as 
in the differential equation derivation we assume a solution of the 
form 


w(x,t) - ^2 ^ n (x) (A. 1-25) 

n=l 

Finally, the equations of forced motion are readily simplified to 

2 C % 

MjUt) +MjWj5j(t) = / F z (x,t) (x) dx (A. 2-20) 

-i 

where 


M . 
: 



j (x)m(x) 


dx . 


(A. 2-21) 


One mode relates the rigid body translation and another the rigid body 
rotation. Letting these be the first and second modes, their corres- 
ponding mode shapes are 


4*1 (x) “ 1 (A. 2-21) 

4> 2 ( x ) = x- (A. 2-22) 

The derivation of the analogue of Eq. (A. 2-20) is given for the two 
dimensional case in Section 5. There is no difference in approach 
between the one and two dimensional situations, hence the proof of Eq. 
(A. 2-20) is omitted here. 

A* 3 . Free Vibration of the Unrestrained Aircraft 

From Chapter 1, Section A. 2, we have for the aircraft in the 
absence of external forces 
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ff w(x,y,t) p(x,y) dxdy = 0 

(A. 3-1) 

ff w(x,y,t) xp (x,y) dxdy = 0 

(A. 3-2) 

ff w (x,y , t) y p (x, y) dxdy = 0 • 

s 

(A. 3-3) 

In addition, the relation between inertial and elastic 

forces yields 

, 3w (0 , 0 , t) „ 3w(0,0,t) 

w(x,y,t) - w (0 , 0 , t) - x - y 


= ~ ff c(x,y;C,n) p(C,n)w (C,n,t) d£dn. 

(A. 3-4) 

Upon imposition of a solution of the form 


w(x,y,t) = W(x,y)T(t) 

(A. 3-5) 

the four equations of motion become 


ff W(x,y)p(x,y) dxdy = 0 

s 

(A. 3-6) 

ff W(x,y) xp (x,y) dxdy = 0 

(A. 3-7) 

ff W (x,y) yp(x,y) dxdy = 0 

s 

(A. 3. 8) 

T [W(x,y) -W(O.O) 1 


= ~^ff C(x,y; C»h) p(£,u) W(C,n) dCdn. 

(A. 3-9) 

Note that Eq. (A. 3-9) is separable and reduces to 


T +uj 2 T = 0 

(A. 3-10) 
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W (x, y) - W(0,0) -x u - -^' 0) 


.. 3W(0, 0) 

Y 3y 


= CJ 2 J C-(x,y; ?,n>p(£/n)W(£,ri)d£dri (A. 3-11) 

S 

where aj 2 is a constant of proportionality of dimension sec -2 . 

As was the case in the one dimensional beam theory, we shall be 
able to eliminate the terms representing the amplitude and deflection 
at the center of gravity. Multiplying Eq. (A. 3-11) by p(x,y) and 
then integrating over the surface of the airplane, we obtain 


If W(x,y)p(x,y) dxdy - W(0,0) If p(x,y) dxdy 


_3 W ( 
3x 


xp(x f y) dxdy yp(x,y) dxdy 



C(x,y;5,n) 


P(?fh) W(C,n) dCdndxdy. 


(A. 3. 12) 


But the first integral vanishes by Eg. (A. 3-6) and the third and fourth 
integrals likewise vanish by the definition of center of gravity. 

Hence Eg. (A. 3-12) may be solved for W(0,0) as follows; 


M SS s p c ^ x,y; ^ ,r l) P (£,h) W(£, n) d£dndxdy 

(A. 3-13) 

where M is the mass of the aircraft. 

Returning to Eg. (A. 3-11), multiplying by xp(x,y) followed by 
integration over the surface yields 


3W(0 ,Q) 
3x 


“ " T-ff^9(x,y) Jf C(x,y;C,n 
y s 


) p(C,n) 


where 


I 


y 



(x,y)dxdy. 


w ( £/ n) d£dr]dxdy 
(A. 3-14) 
(A. 3-15) 
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Similarly, multiplication of Eq. (A. 3-11) by yp(x,y) and subsequent 
integration leads to 

3W( P' 0) - = ff yp(x,y) > /jTc(x,y; 5 ,n) p(€/n)W(C,n)dCdndxdy 

y x s s 

(A. 3-16) 


where 


I 


x 



2 


y p 


(x,y) 


dxdy . 


(A. 3-17) 


Substituting Eqs . (A. 3-13), (A. 3-14) and (A. 3-16) into Eq . (A. 3-11) 

we obtain the integral equation governing the mode shapes of the free 
vibration of the aircraft 


W(x,y) 



G(x,y;£,n)W(C,n) 


p(C,n)d£dn 


(A. 3-18) 


where 


G(x,y;5,h) =C(x,y;£,n) -JJ C (r , s ; £, n )J^ p(r,s)drds. 


(A. 3-19) 


A. 4 . Orthogonality Condition for Mode Shapes 

The mode shapes of the unrestrained aircraft can be shown 
orthogonal as follows. For any two different modes we can write, from 
Eq . (A. 3-18) 


W n (x,y) 

w m (x,y) 


w 3 /T G(x,y; C,h)W n (Crn) P (Cf h) <3£dr| 

G(x,y;£,n) W m (£, n) P U, n) d£,dn . 
s 


(A. 4-1) 


(A. 4-2) 
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Multiplying Eqs. (A. 4-1) and (A. 4-2) by W m (x,y) p (x,y) and Yl^ (x,y) p (x,y) 
respectively and integrating over the surface of the airplane, we 
obtain 


X£T W n W ro PdXdy = ff s "™ p [ff s GW n pd5dri dxdY 


(A. 4-3) 


hfj V 

ms S 



GW m pd^dr)dxdy . 


(A. 4-4) 


By interchanging the dummy arguments and then reversing the order of 
integration, the right side of Eg. (A. 4-4) becomes 


pdxdy 


~ I I W W 

a ?JJ m n 

m s 

fj W m (x ' Y) p(x ' y) [yy"G(C,n;x,y) W n (£, n) P (C ,n) dCdnjdxdy. 


(A. 4-5) 


In Section 3 the influence function for the unrestrained 
airplane was derived 


G(x,y;£,n) C(x,y;C,b) C(r,s;g,n)[~ ^+ Xr + y s l p (r , s) drds . 

J *s L y x J 

(A. 3-17) 

Noting that C(x,y;£,n) = C ( £ ,n ; x ,y),[See Ref.l, Chapter 2] it is 
obvious from Eg. (A. 3-17) that 


G(S,n?x,y) = G(x,y;C,n) + ffc (r , s ; C , n) p (r , s) drds . 

* •'g *- A y i x J 

~ ff C(r,s;x,y) Jif|£+^£]p( r ,s)drds. 


(A, 4-6) 


y x- 
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Substituting Eq. (A. 4-6) into Eq. (A. 4-5) we arrive at 


^ 2 fl ( W m W n pdxdy =// w m (x 'y> P(x,y)|^"[G(x,y;C, 


n) 


+[ f c(r,s;c,n)f^-I £+ F-]p (r ' s) drds 

J J ^ L y x 


-f f C(r f s;x,y)[i+|^+^-]p (r,s) drds] 
J J* y x 


•W n (£, n) P (£, n) d£dn dxdy. 


(A. 4-7) 


Upon expanding Eq. (A. 4-7) into the sum of three intearals , it is 
immediately clear that the third integral vanishes by application of 
Eqs . (A. 3-1), (A. 3-2) and (A. 3-3), as both £ and n appear to the zero 
or first power in the inner integral. If we interchange the order of 
integration in the second integral, performing the integration with 
respect to x and y first, then a similar situation occurs, and by Eqs. 
(A. 3-1), (A. 3-2), and (A. 3-3) again, the second integral will contri- 

bute nothing. Whence Eq. (A. 4-7) simplifies to 


^//" W m W n Pdxdy = f f W m [ GW n P d £ dr ] dxdy . (A. 4-8) 


Comparison with Eq. (A. 4-3) immediately reveals 


( 


1 1 

~Z 2 

oi ai 



w w 


m n 


pdxdy = 0 


(A. 4-9) 


which for m / n is the desired orthogonality condition. 
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Furthermore, the mode shapes corresponding to the rigid transla- 
tion and rotation of the unrestrained aircraft under external forces, 
given by Eqs. (1-15) and (1-16) respectively, are also orthogonal to 
all the deformation mode shapes and each other. This is an immediate 
consequence of Eqs. (1-1) and (1-2) and the definition of the center 
of gravity. 

A. 5 Forced Motion of the Unrestrained Aircraft 

From Chapter 1, by equilibrium of forces and moments, we have 


/ j w(x,y,t) p (x,y) dxdy = JJ F z (x,y,t) dxdy 


(A. 5-1) 


Jj w(x,y,t) xp (x,y) dxdy = JJ F g (x,y, t) xdxdy 


(A. 5-2) 


and relating inertial and elastic forces we obtain 


w(x,y,t) - w (0 , 0 , t) - x- " (0 ' 0,t) 

dX 



C(x,y;r,s) 



(r,s, t) 


w(r,s,t) p (r,s)j drds . 

(A. 5-3) 


Assuming a solution of the form analogous to that for motion of the 
one dimensional beam, 


W(x,y,t) 4> (x,y) £ (t). 

n=l n 


(A. 5-4) 


The rigid body motion can be characterized by the modes 


< ! ) 1 (x,y) =1 oj 1 = 0 


(A. 5-5) 


<f> 2 (x,y) = x w 2 = o. 


(A. 5-6) 
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The other modes represent the elastic deformation with respect to the 
inertial plane of the aircraft, and hence from Eq. (A. 3-11) we have 

34 ) ( 0 , 0 ) 

4> n U,y) - 4> n (0,0) - x ^ 


= c (x,v;r,s) p (r,s) 4> n <r,s) drds (A. 5-7) 

s 

with the additional constraints from Eqs. (A. 3-7) and (A. 3-8) 

Jf <j> (x,y)p(x,y) dxdy = 0 (A. 5-8) 

s 

fj 4> n (x,y) xp (x,y) dxdy = 0 • (A. 5-9) 

s 

The substitution of Eq. (A. 5-4) into Eq. (A. 5-1) in conjunction 
with Eqs. (A. 5-5) and (A. 5-8) yields simply 


£ JJ p(x,y) dxdy = =JJ F z (x,y,t) dxdy. (A. 5-10) 

Similarly, substitution of Eq. (A. 5-4) into Eq. (A. 5-2) together with 
Eqs. (A. 5-6) and (A. 5-9) gives 

I ff x 2 p(x,y) dxdy = M^ 2 =ff F z (x,y , t) xdxdy. (A. 5-11) 

s s 

Note that and M 2 are the mass and pitching moment of inertia 
respectively . 

If we now introduce Eq. (A. 5-4) into Eq. (A. 5-3) we obtain 
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*n ( 0 ' 0 ) 


3 4 >_ ( 0 , 0 ) , 

X - ^x >n (t) 


00 

= ff c < x /y? r /S) [F z (r,s,t) - p(r,s) ^<J> n (r,s) C n (t)l drds. 
s n=l -* 


(A. 5-12) 


Noting that 

3 4> n (0, 0) 

4>j(x,y) - <t>j(0,0) - x ^ = 0 (A. 5-13) 

f° r j = l an d 2, and applying Eq . (A. 5-7) we can rewrite Eq. (A. 5-12) 

9<j> (0,0) 

V x -y> - V°' 0) - x 53E 

”n 

=JJ C(x ' y;r,s) [ F z (r,s,t) - p(r,s)(C 1 (t) + re 2 (t))]drds. 

S 

(A. 5-14) 

Multiplying by <(> m (x,y) p(x,y) and integrating over the surface yields, 
after the application of Eqs . (A. 5-8) and (A. 5-9) and the orthogonality 
condition of the mode shapes, 


E 


e n (t) 

(C n (t) + -a-2-) 


M m + O = ff ^ m < x 'y)p( x »y) ff C(x,y;r,s) • 
m s s 

^F z ( r ,s,t) - p (r,s) + r£ 2 ) J drds. (A. 5-15) 

As the influence function is symmetric, [See Ref. 1, Chapter 2] 

C(x,y;r,s) = C ( r , s ; x , y ) (A. 5-16) 
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we can interchange the variables of integration and the order of 
integration in Equation (A. 5-15) to achieve 


M 

m 


+ ^m) = ff s [V x ' y 't> - p < x ' y) <«i + x ^2 } ] 

m 

• C (x,y; r, s) <|> m (r ,s) p (r, s) drdsdxdy. (A. 5-17) 


The use of Eq. (A. 5-7) allows us to simplify the right side of Eq. 
(A. 5-17) even further 


“»(% + 0 = 4 //h< x ^> - p(x ' y) { 'h + x ^2 ) ]‘ 

'll) ' li)„ J •'s L 

m m 


3<i> (0,0) 

4> m (x,y) - 4> m (0,0) - x Jdxdy. (A. 5-18) 


Expanding the right side of above and making use of Eqs. (A. 5-8), 

I 

(A. 5-9), (A. 5-10) and (A. 5-11), together with the definition of 

center of mass, Eq. (A. 5-18) reduces considerably to 


M m^m + u m^m = /( F z {x '^ t) V X ' y) dxdy ’ (A. 5-19) 

Although the derivation of Eq. (A. 5-19) was restricted to the elastic 
deformation modes, Eqs. (A. 5-10) and (A, 5-11) show that the rigid 
body motion can be similarly represented. 
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APPENDIX B 

Rigid Body Motion - Longitudinal Group 

For this analysis, we will assume the aircraft is in an equili- 
brium constant altitude flight condition. We choose the coordinate 
system as follows: the x-axis coincides with the equilibrium direction 

of flight, the z-axis is chosen down, and the y-axis is taken out the 
wing to form a right handed coordinate system. As usual, the origin 
is at the airplane center of gravity. Once set, the axes are fixed 
to the airplane, and rotate along with the aircraft. 


aircraft 



Figure B1 . Airplane in Equilibrium Flight Condition 


Following Figure Bl, we will denote the aircraft forward velocity 
by 


U = U Q + u (B-l) 

where u represents the velocity perturbation along the x-axis. The 
velocity along the z-axis is given by 
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w 


(B-2) 


where a is the perturbation angle of attack. Note that a positive 
increment in without a rotation of the aircraft induces a positive 
z-component of velocity. We let 0 represent the rotation in inertial 
space of the airframe about the y-axis. Thus 0, the pitch angle, 
measures the rotation of the x-axis from its initial alignment. 
Finally, the deviation from nominal altitude is given by h. 

With the mass and pitching moment of inertia designated by m and 
I respectively, we can write the linearized inertial equations of 
motion 


6X 

6Z 

6P 


d ,U , 

mu 0 dt ( U 0 5 


.. ,da d0. 
mU 0 ( dt 


d 2 0 


dt 


(B-3) 
(B — 4 ) 
(B-5) 


where X and Z denote the forces along the x and z-axes, and P is the 
pitching moment. The first order perturbations in gravity forces 
due to a rotation of the aircraft are given by 


<$X g = -mg 


(B-6) 


<5Z 

g 


= o. 


(B-3) 


Following the usual convention, the aerodynamic forces and moments are 
normalized as below: 
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x a ' 1 “ U 0 SC x <B_8) 

z a * I 4 C z <B ' 9> 

P a - ipU B ScC m (B-10) 


where p is the atmospheric density, S is the aircraft lifting surface 

1 2 . 

area, and c is the reference wing chord length. The term ^ P u 0 1S 
called the dynamic pressure, and will henceforth be designated by q. 

Two other fundamental parameters for any aircraft analysis are 
the coefficients of lift and drag, and C^. The lift coefficient 
is readily calculated by 


C 


L 



(B-ll) 


For supersonic flight, the following approximation is quite accurate 


C 


L 


4a, 


(M 


1 ) 


(B-12 ) 


where M is the Mach number , and allows for a calculation of the 
1 equilibrium angle of attack. Again, for supersonic flight, we can 
approximate the drag coefficient by 


= C, 


1 

? 2 

+ C T 1 2 


3C D n 2 
C_ + 5 — c. 

D 0 d(C^) 


( B-13) 


where C is the skin friction component of drag. 

B r, 
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The normalized forces and moment, obviously zero at equilibrium, 
are functions of the flight variables, U Q , a^, h Q , pitch angle, and 
elevator deflection angle. For small perturbations in these 
quantities, we can write, to first order. 


= C + C a 

x U„ x„ 

u 0 a 


(B-14) 


C 

z 




C 


z 


a + 
a 




h_ 

h„ 


(B-15) 


= C 


m 


u 


2U, 


C C! + a 
m* m 
a a 


c 

2U 


C 0 

o m e 


+ C <5 
m 6 


(B-16) 


where 6 is the perturbation elevator deflection angle. The coeffici- 
ents on the right side above are stability derivatives, partial deri- 
vatives of the forces and moments with respect to equation variables. 

The term C is the speed damping derivative as it relates the 
u 

resistance due to a change in forward speed. For jet aircraft, we 
can use the approximation 



9C x 


3C r 


< 2C D + 


3P 


M) 


(B-17) 


Similarly, the term C z reflects the lift due to an increase in speed 

u 


3C 


8C 

_(2c l + if 1 M) 


The derivative C m is defined as 


ac 


m 


m 




3C 

ii 

3P 


M 


(B-18) 


(B-19) 
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and relates changes in pitching moment to an increase of speed. It 

is primarily due to the center of pressure moving aftward with an 

increase of speed in the transonic region. For high Mach number 

3C D 8c l 

flight, this derivative is negligible, as are the and terms. 

The a derivatives relate changes in the resultant forces and 
moment induced by an increase in angle of attack, 




dc 


'd(cj 



(B-20) 


For supersonic flight however. 


dCL 


d(cJ7 


2 2 
1 = (M Z - 1 )* 



(B-21) 


so that Eq. (B-21) simplifies to 


C = -C T . (B-22) 

X 1j 

a 

The C 

z 

a 


The stability derivative C m relates changes in pitching moment 

a 

due to the change in aerodynamic center of pressure effected by a per- 
turbation in angle of attack. This term is negative for a stable 
aircraft, lest the effect of a small change in angle of attack is to 
induce a motion driving the airplane further from its equilibrium state, 


term is given by 


3C z 

3a 


k 3a 


+ V 


(B-23) 


105 



(B-24) 


m 


3C 3C 3C 

m _ m L 

3a “ 3C 3a ' 

L 


The value of decreases in magnitude as the center of gravity 


moves aft, and for the B-2707 prototype is approximately -0.2. 

The derivative C m . is called the damping in pitch, and represents 


0 

the aerodynamic effects that accompany a rotation of the airplane in 
the xz-plane; 



3C 


m 


3 < 4 § > 


(B-25) 


This term is primarily due to the effect of the tail, and for the SST 
model is approximately -3. 

The a derivative C m . arises since given a change in a, the pres- 

a 

sure distribution over the lifting surfaces does not adjust instan- 
taneously to an equilibrium value. This term is very difficult to 
determine for the wings and body of the airplane as it involves un- 
steady flow. The effect of the tail is derived primarily from the 
concept of the lag of the downwash. The downwash angle, e, measures 
the downward deflection of the airflow at the tail due to the influence 
of the wings. Letting 1 T denote the distance from the center of gravity 
to the center of the tail, and S T the surface area of the tail, we 
can write 



3C 


3 



2 , 1 T.2 S T , 3C L . 3e 
“■ c ~s <Ta T Ja‘ 


(B-26) 


(B-27) 
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In the swept-wing configuration of the B-2702, the C m< term is 

a 

negligible . 

The derivative C denotes changes in the force along the z-axis 
Z P 

due to an increase in altitude, and is given by 



3C_ 


3(^) 


-C 


L* 


(B-28) 


The 6 derivatives relate changes in forces and moments with 
respect to elevator deflections. In the swept-wing configuration, 
and elevator deflection (measured positive downward) has the effect 
of changing the angle of attack for a portion of the lifting surface. 
This creates a force equal to that supplied by a lifting surface of 
area the same as the elevator. Thus 


C 


z 


6 


8C z 



(B-29) 


where S_ is the area of the elevator. Similarly, 
£ 


c 

m 5 


3C 

m 

“37 


Ve 

cs 



(B— 30) 


.where 1^ is the distance from the center of gravity to the elevator. 

As the aircraft is unrestrained, the sum of the aerodynamic and 
gravitational forces must be in equilibrium with the inertial accele- 
ration. From Eqs . (B-3) through (B-10) and (B-14) through (B-16) 

we obtain the equations of motion: 


mU r 


qS U r 


- C 


u 

Or 


a + 


mg 

qS 


0 = 0 


(B-31) 
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(B-32) 



The four equations above are called the airplane longitudinal 
group and are most often presented in the form below, where D stands 
for the differentiation operator with respect to time: 



(B-35) 
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Eq. (B-35) can be put in the form of a minimum realization linear 
dynamical system as follows: 


*1 - * 1±1 T “l u l 


(B-36) 


where 


= [ u a 0 6 h ] 


(B-37) 


Uy m 

U[) q s 


c / — 

x ' qS 
a 


C z 

u ^ 


c */ q s 


C C 

-rr ( -u^2u^m. V q s 


„ ulu n 

q Sc m 4. . c r r / — 


I ' m 


^ n 111 « 

a 0 a a 


2p 

h 0 7 q s 


q qSc c .y. +( -, \ 

X 2U n l nil ^ m* 

y 0 0 a 


„ z mU n 

qSc,_c_ r p/ 0, 

I 2Ut" m- h/qS 
y 0 a 0 
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and the control 


u 


1 


6 . 


(B-40) 
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APPENDIX C 


Derivation of Feedback Control Gains 

Consider the following time invariant linear system with quad- 
ratic cost criterion 


J 


x(t) = Ax ( t) 


x' (T)Qx(T) + 


r 


+ Bu (t) 


[x- (t) u 1 (t) ] 


L 

N’ 


N 

R 



*x(t)‘ 


U (t) 


(C-l) 

(C-2) 


We will show that the control u*{t), which minimized J subject to 
Eq. (C-l), exists, and that this control is a linear functional of 
the state. We shall restrict our model by assuming the following 
conditions on the matrices L, N, R, and Q: 


L' = L ^ 0 (non-negative definite) (C-3a) 

q' = Q 5^ 0 (C-3b) 

R' = R > 0 (positive definite) (C-3c) 

L - NR -1 N' > 0. (C-3d) 


The symmetry and definiteness relations of Eq. (C— 3 ) involve no real 
loss of generality as far as practical, real-world systems go. Any 
quadratic form x'Sx can be put in the form 

x'Sx = x'S^x where S^' = S^. (C-4) 

The definiteness relations guarantee that a minimum for J does exist 
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for all possible state trajectories and controls. Obviously, that 
minimum is zero . 

First, we state and prove a simple lemma. For x and u related 
as in Eq. (C-l) , and K(t) a matrix differentiable on 0,T, then 


T 

K(t)+K(t)A+A'K(t) K (t) B 


x (t) 

[x- (t) u 1 (t) ] 




0 

B'K(t) 0 


u ( t) 


- x' (t) K (t) x (t) 


T 

0 


0 . 


( C— 5 ) 


The proof of Eq. ( D 5 ) is as follows: for any differentiable x and 

K we have 


/ 


[x' (t)K(t)x(t) ] dt - x' (t)K(t)x(t) 


= 0 . 


(C-6) 


Eq. (C-6) can be rewritten 

/V<t)K(t)x(t) + x' (t)K(t)x(T) + x ' (t) K ( t) x ( t) ] dt 


- x' (t)K(t)x(t) 


= 0 . 


(C-7) 


Substitution of Ax(t) + Bu(t) for x(t) leads immediately to Eq.(C-5). 

Now we will show that the control which minimizes J is given by 


u*(t) = -R -1 (N' + B'K(t))x(t) 


(C-8) 


where K(t) satisfies the differential equation 


K(t) = -K(t) A-A'K(t) + (K(t)B+N)R -1 (N'+B , K(t) ) -L 


(C-9) 
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with terminal condition K(T) = Q * . It is clear from Eq. (C-9) that 
the matrix K(t) , if it exists, is symmetric. The definiteness 
restrictions of Eq. (C-3) insure that a positive definite solution of 
Eq. (C-9) will exist. 

If we add the identity of Eq. (C-5) to the cost functional of 
Eq. (C-2) , we obtain 


J 



- 


“ ~ 

L N 


X 

N' R 


u 

_ _ 


- . 


dt + x' (T)Qx(t) 



u'] 


K+KA+A'K 


B'K 



x'Kx 


T 

0 


(C-10) 


Combining the two integrals and using the differential equation for 
K(t) yields 




(KB+N) R _1 (N'+B'K) 

N+KB 

X 

P [X' 
0 

U'] 

N'+B 'K 

R 

u 





„ _ 


+ x' (0) K (0) x (0) . 


(C-ll) 


The integral term in Eq. (C-ll) is obviously non-negative as 
follows from Eq. (C-3) . However if u(t) is chosen according to Eq. 
(C-8) , then the integral term is identically zero and the total cost 
function simplifies to 

*See note at end of Appendix. 
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J 


(C— 12 ) 


= x(0)K{0)x{0) . 

Not only does the feedhack gain K(t) figure in the calculation of the 
control, but it gives an explicit formulation for the cost of the 
optimal trajectory. 

With the minimum cost control defined by Eq. (C-8) , the 
optimal state trajectory is given by 

x(t) = A -BR" 1 (N 1 + B'K(t)) x(t). (C-13) 


Note that Eq. (C-9) can be rewritten as 


K = -K(A-BR V) - (A-BR -1 N' ) 'K +KBR~ 1 B , K - (L-NR _1 N ' ) , 

(C-14) 

and it is easily shown that the problem defined by Eqs. (C-l) and 
(C-2) is equivalent to the problem 


min J = 


x' (T)Qx(T) 



u’] 


L-NR ^N ' 


0 


x 


dt 


R 


u 


subject to 


(C-15) 


x(t) = (A - BR _1 N ' ) x (t) + Bu(t) . 


(C-16) 


This last assertion can be established by an argument paralleling that 
used to derive the minimum cost control of the original system. 


*[Note: The reader, should note that there has been no justifi- 

cation for the boundary condition K(T)=Q. This is correct, however, 
and can be deduced by other means TnclucTing transversality conditions.] 
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APPENDIX D 


Solution of the Algebraic Riccati Equation 
Consider the linear system and quadratic cost functional 

x = Ax + Bu; x(0) = x Q (D-l) 

dt. (D-2) 

We have already seen that the minimum cost control is given by 

u* = -R~ 1 B 1 Kx (D-3) 

where K is the solution of the matrix quadratic equation 

A'K + KA -KBR -1 B'K + C = (D-4) 

For a simplification of notation, let us define 

D = BR -1 B'. (D-5) 



If x is an n-dimensional vector, then the matrices A, C, D, 
and K will be n x n dimension matrices . We will show that the 
solution to Eq. (D-4) can be expressed in terms of the eigenvectors 
of the 2n x 2n matrix 


V 




(D-6) 
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Further, the poles of the closed loop minimum cost trajectory, 

x = (A - BR -1 B'K) x (D-7) 

coincide with the eigenvalues of V. 

It will be assumed that the matrix V has a diagonal Jordan 
canonical form. This will be true if the eigenvalues of V are distinct, 
or in the case that this fails, there exist 2n linearly independent 
eigenvectors. The method to be employed can be extended to cases 
where V is non-diagonalizable. 

We will first establish that the eigenvalues of V have 
quadrangular symmetry; i.e. if X is an eigenvalue, then so are X, -1, 
and -X (where the superbar denotes complex conjugate) . For the 
matrices A, C, and D, defined over the field of real numbers, it is 
clear that complex eigenvalues will always occur in conjugate pairs. 
Consider the 2n x 2n matrix 


J 



I 



(D-8) 


Notice that J 


2 


-I. Thence, as is verified by direct multiplication 


JVJ = V' 


(D-9) 


which leads immediately to 
-JV = V'J. 


(D-10) 
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If X is an eigenvalue of V with associate eigenvector e, then from 
Eq. (D-10 ) 

V'Je = -JVe = -XJe (D-ll) 


which shows that -X is an eigenvalue of V* with associated eigenvector 
Je. But as V and V' have the same eigenvalues, we have shown that -X 
is an eigenvalue of V. It can be further shown that for a well-posed 
system 

x = Ax + Bu ; y = F ' x (D-12) 


where FF'=C, completely controllable and completely observable, no 
eigenvalue of V will be purely imaginary. (It should be noted that 
for the case of real valued eigenvalues, the quadrangular symmetry 
condition collapses to a simple symmetry condition; if X is a real 
eigenvalue, then -X is also an eigenvalue.) 

Thus we can take the n eigenvalues with positive real parts, 
Xl'..., x n ( the other eigenvalues are -X^..., -X n ) and arrange their 
corresponding eigenvectors in column form [e^, ..., £ n ] . This 2n x n 
matrix can be partitioned into two n x n matrices X and Y by 


X 

Y 


[e L , 


e ] 
— n 


(D-13) 


Hence we can write 




- 


— 

A' C 


X 


X 

D -A 


Y 


Y 



— 



- — 


- 


- - 


(D-14) 
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where A is the n x n Jordan form for the eigenvalues of V with 
positive real parts. The solution to Eq. (D-4) may be written 


K • = 


XY 


(D-15) 


The proof is quite straightforward. From Eq. (D-14) we have 


A'X + CY = XA (D-16) 

DX - AY = YA . (D-17) 

Postmultiplying Eqs . (D-16) and (D-17) by Y -1 (it can be shown that 

Y is invertible) , and premultiplying Eq. (D-17) by XY -1 yields 

A'XY -1 + C = XAY -1 
XY~ 1 DXY~ 1 - XY -1 A = XAY -1 . 


(D-18) 

(D-19) 


Subtraction of Eq. (D-19) from Eq. (D-18) results in 


A' (XY -1 ) 


(XY 1 ) A 


, -1 

(XY ) D ( XY ■*■ ) + C = 


0 . 


(D-20) 


It can be further established that for C and D symmetric and positive 
semi-definite, with at least one of C or D non-singular, then K will 
be positive definite. 

To show that the eigenvalues of V are the poles of the closed 
loop system, let us define the matrix G by 



(D-21) 


Then post multiplication of Eq. (D-17) by Y -1 yields 
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DK 


A 


G. 


(D-22) 


Since Eq. (D-21) represents a similarity transformation, the eigenvalues 
of G and A are identical; hence the eigenvalues of A - DK, which are 
the poles of the closed loop minimum cost trajectory by Eq. (D-7) , are 
given by the diagonal elements of -A. Furthermore, if all the eigen- 
values comprising A have positive real parts, then the system of Eq. 

(D-7) will be asymptotically stable. 
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